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Figure  1.  The  distribution  of  f-k  azimuth  as  a  function  of  signal-to-noise  ra¬ 
tio  is  plotted  for  the  Pn,  Pg,  Sn,  and  Lg  phases  of  events  from  the  mining 
area  north  of  Apatity  (Kola  Peninsula),  located  at  ARCESS.  For  each  plot, 
the  dotted  line  denotes  the  mean  of  the  observations,  while  the  dashed  line 
corresponds  to  the  actual  azimuth.  The  lines  marked  by  squares  represent  2 
standard  deviations,  as  fitted  by  an  exponential  decay  with  constant  rate . 12 

Figure  2.  The  azimuth  bias  (f-k  azimuth  minus  true  azimuth)  is  plotted  for 
the  four  regional  phases,  Pn,  Pg,  Sn,  Lg,  from  events  corresponding  to  the 
mine  locations  in  Table  1  and  recorded  at  ARCESS.  The  circles  denote  neg¬ 
ative  bias,  and  the  crosses  denote  positive  bias.  The  symbol  size  is  propor¬ 
tional  to  the  amount  of  bias.  The  location  of  the  ARCESS  array  is  shown  by 
a  triangle.  The  arc  represents  a  distance  of  500  km  from  the  array . 14 

Figure  3.  The  standard  deviation  of  the  f-k  azimuth  for  a  signal-to-noise  ra¬ 
tio  of  10  is  plotted  for  the  four  regional  phases,  Pn,  Pg,  Sn,  Lg,  from  events 
corresponding  to  the  mine  locations  in  Table  1  and  recorded  at  ARCESS. 

The  symbol  size  is  proportional  to  the  standard  deviation  in  (the  maximum 
is  10 ).  The  location  of  the  ARCESS  array  is  shown  by  a  triangle.  The  arc 
represents  a  distance  of  500  km  from  the  array . 16 

Figure  4.  The  observations  of  azimuth  bias  in  Figure  2  were  interpolated  and 
extrapolated  to  cover  a  regular  latitude  and  longitude  grid  (with  increments 
0.5  in  longitude  and  0.25  in  latitude).  The  resulting  values  (in  )  are  displayed 
on  contour  plots  for  the  four  regional  phases,  Pn,  Pg,  Sn,  and  Lg . 17 

Figure  5.  The  observations  of  azimuth  scatter  in  Figure  3  (standard  devia¬ 
tion  in  for  an  SNR  of  10)  were  interpolated  and  extrapolated  to  cover  a  reg¬ 
ular  latitude  and  longitude  grid  (with  increments  0.5  in  longitude  and  0.25 
in  latitude).  The  resulting  values  are  displayed  on  contour  plots  for  the  four 
regional  phases,  Pn,  Pg,  Sn,  and  Lg . 18 

Figure  6.  The  azimuth  bias  is  plotted  for  Pn  from  events  corresponding  to 
the  mine  locations  in  Table  1  and  recorded  at  ARCESS.  At  the  top,  the  dif¬ 
ference  between  the  f-k  azimuth  and  the  mine  azimuth  is  plotted,  while  at 
the  bottom  is  shown  the  difference  between  the  azimuth  from  polarization 
of  3-component  data  (from  Suteau-Henson,  1992b)  and  the  mine  azimuth. 

The  circles  denote  negative  bias,  and  the  crosses  denote  positive  bias.  The 
symbol  size  is  proportional  to  the  amount  of  bias.  The  location  of  the 
ARCESS  array  is  shown  by  a  triangle.  The  arc  represents  a  distance  of  500 
km  from  the  array . 19 
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Figure  7.  The  phase  identification  performance  (percentage  of  correctly  identified 
detections)  is  displayed  for  eleven  polarization  recipes  (described  in  the  text  and 
Table  3)  at  each  of  five  IMS/IRIS  stations.  Detections  in  the  three  groups  regional 
S,  regional  P,  and  teleseismic  were  classified.  For  each  station,  the  standard  recipe 
(S),  the  optimum  recipe  for  azimuth  estimation  from  Riviere-Barbier  et  al.  (1992), 
when  it  applies,  and  the  preferred  recipe  from  this  experiment  are  highlighted . 27 

Figure  8.  The  phase  identification  performance  (percentage  of  correcdy  identified 
detections)  is  displayed  for  eleven  polarization  recipes  (described  in  the  text  and 
Table  3)  at  station  ARU.  Detections  in  the  three  phase  groups,  regional  S,  regional 
P,  and  Teleseismic,  were  classified.  In  spite  of  somewhat  lower  overall  perfor¬ 
mance  (see  Figure  7),  recipe  A3  is  preferred  to  K3  due  to  its  higher  success  rate  for 
the  regional  P  group . 28 

Figure  9.  The  phase  identification  performance  (percentage  of  correctly  identified 
detections)  is  displayed  for  each  of  the  four  groups  (regional  S,  regional  P,  Teleseis¬ 
mic,  and  Noise),  and  their  average,  at  each  of  five  IMS/IRIS  stations.  The  results 
are  shown  for  the  standard  recipe  “S”,  the  optimum  recipe  “A”  for  azimuth  estima¬ 
tion  (from  Riviere-Barbier  et  al.,  1992)  when  it  applies,  and  the  best  recipe  “B”  ob¬ 
tained  for  the  S  /  P  /  T  classification  problem . 30 

Figure  10.  The  distributions  of  azimuth  residuals  (azimuth  from  polarization  minus 
azimuth  from  NEIC  location)  of  teleseismic  P  arrivals  are  displayed  for  eleven  po¬ 
larization  recipes  (described  in  the  text  and  Table  3),  at  each  of  four  IMS/IRIS  sta¬ 
tions . 31 

Figure  11.  The  combined  performance  for  phase  identification  and  azimuth  estima¬ 
tion  is  displayed  for  a  set  of  eleven  polarization  recipes  (described  in  the  text  and 
Table  3),  at  four  IMS/IRIS  stations.  The  classification  problem  considered  is  the 
separation  of  regional  S,  regional  P,  and  teleseismic  detections . 33 

Figure  12.  The  combined  performance  for  phase  identification  and  azimuth  estima¬ 
tion  is  displayed  for  four  of  the  IMS/IRIS  stations.  The  data  set  used  for  phase  iden¬ 
tification  included  all  four  detection  groups  (regional  S,  regional  P,  Teleseismic, 
and  Noise).  At  each  station,  the  performance  was  estimated  for  the  standard  recipe 
S,  the  optimum  recipe  for  azimuth  estimation  from  Riviere-Barbier  et  al.  (1992), 
when  it  applied,  and  a  selected  subset  of  recipes  that  performed  well  for  the  S  /  P  / 

T  combined  optimization  experiment  (see  Figure  11).  The  resulting  preferred  reci¬ 
pes  are  highlighted . 34 
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Figure  13.  The  map  shows  the  geographical  distribution  of  the  events  and 
the  seismic  stations  for  which  waveform  data  were  analyzed.  The  two  un¬ 
derground  nuclear  explosions  in  the  South  Pacific  are  marked  with  cross¬ 
es . 38 

Figure  14.  The  signal  processing  and  transformation  in  the  time  domain,  up¬ 
per  frame,  and  in  the  depth  domain,  lower  frame.  The  time  domain  process¬ 


ing  involves  bandpass  filtering  to  enhance  the  signal-to-noise  ratio,  forming 
of  a  STA/LTA  trace,  masking  of  window  for  expected  arrival  time  of  PcP, 
and  transformation  into  a  binary  trace.  The  binary  trace  in  the  time  domain 
is  mapped  into  two  binary  traces  in  the  depth  domain,  one  each  for  pP  and 
sP . 40 

Figure  15.  The  scatter  diagrams  shows  the  ratios  of  the  periods  of  pP/P,  to 

the  left,  and  of  sP/P  to  the  right,  as  a  function  of  the  delay  times . 41 

Figure  16.  STA/LTA  traces  plotted  as  a  function  of  epicentral  distance  for 
all  the  waveforms.  Only  the  parts  of  the  traces  with  the  STA/LTA  above  the 
threshold,  1.75,  were  marked  as  black  dots . 43 


Figure  17.  The  figure  shows  waveforms  and  binary  traces  sections  for  three 
events  in  the  Afghanistan-USSR  border  region  on :  (a)  91/05/13  at 
23:55:33.7,  (b)  91/04/26  at  22:24:05.3,  and  (c)  91/05/13  08:32:43.8.  Stacion 
codes  and  epicentral  distances  are  indicated  to  the  right  of  each  waveform. 
The  waveform  sections  also  include,  at  the  top,  envelope  stacks  of  pP  and 
sP  as  a  function  of  depth  and  the  standard  deviation  of  the  first  arrival  time 
P  residuals  as  a  function  of  depth  (SIGMA  P  RES  in  the  diagrams).  The  sig¬ 
nal-to-noise  ratio,  SNR,  and  the  P  residual  for  the  NEIS  hypocenter  solution 


are  shown  to  the  right  for  each  waveform . 44 

Figure  18.  (a)  The  figure  shows  how  the  candidate  depth  is  defined  from  a 

pP  depth  trace  stack . 47 

Figure  18.  (b)  The  barplot  shows  the  number  of  events  with  and  without 


candidate  depth  as  a  function  of  depths  reported  by  NEIS.  Most  of  the  NEIS 
depths  were  constrained  to  33  km  and  a  significant  portion  of  those  had  can¬ 


didate  depths . 48 

Figure  19.  The  geographical  distribution  in  the  Western  (a)  and  Eastern  (b) 
hemisphseres  of  events  with  and  without  candidate  depths . 49 
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Figure  20.  Success  rates,  given  as  percentages,  to  define  candidate  depths  as  a  func¬ 
tion  of  the  number  of  waveforms  and  azimuthal  coverage.  Azimuthal  coverage  is 
given  as  the  number  of  45  degree  sectors  around  the  source  to  station  azimuth  for 
which  there  is  one  or  more  station  with  waveform  data. . 51 

Figure  21.  Depths  taken  as  the  average  of  the  lower  and  upper  values  of  the  can¬ 
didate  depths  (‘pP  DEPTHS’  in  the  diagrams)  compared  with  NEIS  depths.  Com¬ 
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Mapping  Azimuth  Bias  and  Scatter  Around  ARCESS 


Anne  Suteau-Henson 

Introduction 

Azimuth  biases  of  events  located  by  arrays  are  well  documented  (e.g.,  Mykkeltveit 
et  al. .,  1990;  Vogfiord  and  Langston,  1992).  To  our  knowledge,  however,  no  attempt  has 
been  made  to  use  such  knowledge  to  appropriately  correct,  and  thereby  improve,  event 
locations  by  automated  systems  such  as  the  Intelligent  Monitoring  System  (IMS,  Bache  et 
al.,  1990). 

In  a  previous  Report  (Suteau-Henson,  1992b),  variations  of  azimuth  residuals, 
from  both  polarization  and  f-k  analysis,  were  studied  as  a  function  of  azimuth  and  distance 
for  events  recorded  at  the  NORESS  and  ARCESS  arrays  by  the  IMS.  In  another  study 
(Suteau-Henson,  1992a),  the  azimuth  bias  and  scatter  were  systematically  measured  for 
regional  phases  of  events  recorded  at  ARCESS.  The  latter  study  is  extended  here  by  map¬ 
ping  the  observed  bias  and  scatter.  They  apparently  are  slowly  varying  in  the  area  sur¬ 
rounding  the  array,  and  can  be  approximated  by  a  smooth  surface,  through  interpolation 
and  extrapolation  of  die  observations.  The  result  is  tables  of  azimuth  corrections  and  azi¬ 
muth  weights  (inversely  proportional  to  the  amount  of  scatter  for  a  given  signal-to-noise 
ratio)  for  each  regional  phase,  covering  a  latitude  and  longitude  grid  around  ARCESS. 
Such  tables  can  then  be  used  to  improve  automated  locations  by  the  IMS,  in  the  same  way 
traveltime  correction  tables  are  used  to  improve  distance  estimates. 

Observations  of  Azimuth  Bias  and  Scatter 

A  model  was  proposed  by  Suteau-Henson  (1992a),  to  fit  the  statistical  distribution 
of  azimuth  estimates  from  f-k  analysis  for  regional  phases  from  specific  event  clusters 
around  ARCESS.  It  was  shown  that  for  events  associated  with  a  given  mine  (or  mining 
area),  the  distribution  of  azimuth  as  a  function  of  signal-to-noise  ratio  (SNR)  could  be 
described  by  its  mean  (usually  different  from  the  “true"  mine  azimuth  and  phase-specific), 
and  by  its  SNR-dependent  standard  deviation  (usually  phase-specific).  The  SNR-depen- 
dence  of  the  standard  deviation  could  be  fitted  by  an  exponential  decay.  Such  fits  to  the 
data  are  shown,  as  an  example,  for  events  from  the  mining  area  north  of  Apatity  in  Figure 
1.  To  simplify  modelling  of  the  standard  deviation,  a  constant  decay  rate  was  assumed  for 
all  seismic  phases  and  the  entire  area  around  ARCESS.  Such  a  simple  model  fits  the  data 
reasonably  well. 
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Figure  1:  The  distribution  of  f-k  azimuth  as  a  function  of  signai-to-noise  ratio  is 
plotted  for  the  Pn,  Pg,  Sn,  and  Lg  phases  of  events  from  the  mining  area  north  of 
Apatity  (Kola  Peninsula),  located  at  ARCESS.  For  each  plot,  the  dotted  line 
denotes  the  mean  of  the  observations,  while  the  dashed  line  corresponds  to  the 
actual  azimuth.  The  lines  marked  by  squares  represent  ±2  standard  deviations,  as 
fitted  by  an  exponential  decay  with  constant  rate. 
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The  data  set  used  for  this  study  covers  two  years  of  IMS  operation  (from  Novem¬ 
ber  1, 1990  to  November  1, 1992).  The  locations  of  die  mines  (or  mining  areas)  for  which 
azimuth  bias  and  scatter  were  thus  measured  are  given  in  Thble  1 .  The  locations  are  those 


Table  1:  Mine  Locations 


Mine  or  Mining 
Area 

Latitude  N 

Longitude  E 

Kiruna® 

67.833 

20.205 

Malmberget  /  Aitik 

67.15 

20.785 

Saattopora 

67.79 

24.43 

Varangerb 

70.48 

28.50 

Nikel 

69.42 

30.80 

Olenogorsk 

68.16 

33.18 

Apatity 

67.645 

33.87 

Kovdor 

67.56 

30.44 

Elijaervi 

65.78 

24.70 

a.  A.  Dahle  (personal  communication) 

b.  S.  Mykkeltveit  (personal  communication) 


from  the  Helsinki  University  (M.  Uski,  personal  communication),  unless  noted  otherwise. 

Event  clusters  that  were  not  in  the  vicinity  of  a  known  mine  or  quarry  were  not 
included.  Clusters  of  events  that  could  be  associated  with  any  of  several  mines  spanning  a 
wide  azimuth  range  could  not  be  used.  For  mining  areas  with  several  closely  spaced  mines 
(e.g.,  Malmberget/Aitik,  Apatity),  the  average  location  is  given  in  Table  1. 

Mapping  the  Observations 

In  Figure  2  the  observations  of  azimuth  bias  (observed  azimuth  minus  “true”  azi¬ 
muth)  at  the  nine  mine  locations  in  Table  1  are  plotted  on  a  map  of  the  area  around 
ARCESS  for  each  of  the  four  regional  phases,  Pn,  Pg,  Sn,  and  Lg.  Although  the  resulting 
distribution  differs  with  each  phase  type,  some  general  trends  are  observed: 

•  The  azimuth  bias  tends  to  be  negative  southwest  and  south  of  the  array,  and  pos¬ 
itive  east,  northeast  and  southeast  of  the  array; 
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Figure  2:  The  azimuth  bias  (f-k  azimuth  minus  true  azimuth)  is  plotted  for  theiour 
regional  phases,  Pn,  Pg,  Sn,  Lg,  from  events  corresponding  to  the  mine  locations 
in  Table  1  and  recorded  at  ARCESS.  The  circles  denote  negative  bias,  and  the 
crosses  denote  positive  bias.  The  symbol  size  is  proportional  to  the  amount  ol 
bias.  The  location  of  the  ARCESS  array  is  shown  by  a  triangle.  The  arc  represents 
a  distance  of 500  km  from  the  array. 
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•  The  azimuth  bias  varies  smoothly  over  die  area  around  ARCESS,  remaining 
within  10°. 

Similar  maps  are  shown  in  Figure  3  for  the  observed  azimuth  scatter.  As  the  scat¬ 
ter  is  a  function  of  SNR,  the  values  plotted  arbitrarily  correspond  to  the  standard  deviation 
(in  °)  for  an  SNR  of  10. 

Interpolating  the  Observations 

As  the  variations  of  azimuth  bias  and  scatter  around  ARCESS  are  generally  slow, 
they  can  be  fitted  by  a  smooth  surface.  Interpolation  and  extrapolation  of  the  observations 
mapped  in  Figures  2  and  3  were  performed,  resulting  in  values  of  azimuth  bias  and  scatter 
covering  a  regular  latitude  and  longitude  grid  around  ARCESS.  The  contour  plots  of  Fig¬ 
ures  4  and  5  display  these  values  for  the  bias  and  scatter,  respectively.  Note  that  the 
results  for  the  Sn  phase  are  less  reliable,  due  to  sparse  observations.  The  figures  show  that: 

•  The  azimuth  bias  varies  smoothly  from  about  -10°  to  the  southwest  to  about 
+10°  to  the  northeast 

•  Apart  from  a  low  near  Saattopora  (standard  deviation  less  than  3°),  the  azimuth 
scatter  is  characterized  by  largo*  gradients  for  the  P-type  than  for  the  S-type 
phases.  Pn  to  the  north  and  Pg  to  the  northeast  have  large  standard  deviations 
(from  5  to  12°).  The  Lg  scatter  is  larger  to  the  south  (more  than  5°).  Elsewhere, 
the  standard  deviation  is  less  than  4°  for  all  phases  at  this  SNR  level.  Note,  how¬ 
ever,  that  for  a  given  magnitude,  the  SNR  of  secondary  phases  is  always  rela¬ 
tively  low  (as  it  is  measured  against  the  coda),  so  that  on  average,  standard 
deviations  for  S-type  phases  are  typically  largo-  than  those  shown  in  Figure  5. 

Comparison  of  F-K  and  3-C  Pn  Azimuth  Bias 

In  an  attempt  to  obtain  corroborating  evidence  for  the  observed  azimuth  bias  from 
f-k  analysis,  the  azimuth  bias  from  polarization  analysis  of  three-component  data  was  sim¬ 
ilarly  measured  for  the  Pn  phase  at  the  locations  in  Table  1.  The  data  set  used  was  that 
from  Suteau-Henson  (1992b)  which  is  about  four  times  smaller  than  that  used  for  the  f-k 
studies  (hence,  resulting  in  greater  uncertainties)  and  not  overlapping  with  it  The  f-k  and 
3-component  (3-C)  Pn  azimuth  biases  are  compared  in  Figure  6  .  There  is  a  good  qualita¬ 
tive  agreement  between  the  two  sets  of  observations. 
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Figure  3:  The  standard  deviation  of  the  f-k  azimuth  for  a  signai-to-noise  ratio  ot 
10  is  plotted  for  the  four  regional  phases,  Pn,  Pg,  Sn,  Lg,  from  events  correspond¬ 
ing  to  the  mine  locations  in  Table  1  and  recorded  at  ARCESS.  The  symbol  size  is 
proportional  to  the  standard  deviation  in  °  (the  maximum  is  1CP).  The  location  of 
the  ARCESS  array  is  shown  by  a  triangle.  The  arc  represents  a  distance  of  500 
km  from  the  array. 
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Figure  4:  The  observations  of  azimuth  bias  in  Figure  2  were  interpolated  and 
extrapolated  to  cover  a  regular  latitude  and  longitude  grid  (with  increments  0.5*  in 
longitude  and  0.25P  in  latitude).  The  resulting  values  (in  °)  are  displayed  on  con¬ 
tour  plots  for  the  four  regional  phases,  Pn,  Pg,  Sn,  and  Lg. 
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Figure  5:  The  observations  of  azimuth  scatter  in  Figure  3  (standard  deviatioh'in  ° 
for  an  SNR  of  10)  were  interpolated  and  extrapolated  to  cover  a  regular  latitude 
and  longitude  grid  (with  increments  0.5°  in  longitude  and  0.2SP  in  latitude).  The 
resulting  values  are  displayed  on  contour  plots  for  the  four  regional  phases,  Pn, 
Pg,  Sn,  and  Lg. 
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Figure  6:  The  azimuth  bias  is  plotted  for  Pn  from  events  correspornSng  to  the 
mine  locations  in  Table  1  and  recorded  at  ARCESS.  At  the  top,  the  difference 
between  the  f-k  azimuth  and  the  mine  azimuth  is  plotted,  while  at  the  bottom  is 
shown  the  difference  between  the  azimuth  from  polarization  of  3-component  data 
(from  Suteau-Henson,  1992b)  and  the  mine  azimuth.  The  circles  denote  negative 
bias,  and  the  crosses  denote  positive  bias.  The  symbol  size  is  proportioned  to  the 
amount  of  bias.  The  location  of  the  ARCESS  array  is  shown  by  a  triangle.  The  arc 
represents  a  distance  of 500  km  from  the  array. 
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Discussion 

Our  observations  of  azimuth  bias  around  ARCESS  are  in  agreement  with  the  study 
by  Vogfjord  and  Langston  (1992).  They  carefully  analyzed  the  ARCESS  recordings  of  a 
few  large  events  from  known  mines.  They  found  biases  for  the  azimuths  from  f-k  analysis, 
that  generally  depended  on  phase  type.  For  example,  PmP  and  SmS  of  events  from  the 
quarry  on  the  Varanger  peninsula  showed  an  azimuth  5  larger  than  that  of  the  Pg  and  Sg 
crustal  phases.  The  phase  identification  used  in  this  study  is  that  currently  in  use  in  the 
IMS  and  does  not  include  Moho  reflections.  For  close  events  such  as  those  from  the 
Varanger  quarry,  the  f-k  analysis  window  for  the  Pg  phase  includes  both  the  crustal  phase 
(Pg)  and  the  intensive  Moho  reflection  (PmP),  which  explains  why  the  P-phase  azimuth 
has  a  large  positive  bias  (Vogfjord,  personal  communication).  A  similar  explanation  can 
be  given  for  the  S-phase  bias.  Vogfjord  and  Langston  (1992)  speculate  that  the  observed 
biases  are  due  to  lateral  heterogeneities,  that  affect  phases  differently,  depending  on  their 
path  and  frequency  content.  Thurber  et  al.  (1991)  modelled  the  amount  of  bias  produced 
by  different  kinds  of  heterogeneities  (dipping  Moho,  velocity  gradients  at  different  depths 
in  the  crust  and  upper  mantle),  as  well  as  the  distance  range  and  azimuth  span  affected.  A 
better  understanding  of  the  relation  between  laterally  heterogenous  structure  and  azimuth 
bias  would  help  predict  bias  in  areas  with  known  structure,  but  sparse  seismicity  or  little 
“ground  truth”. 

In  Fennoscandia,  more  knowledge  has  been  acquired  about  location  of  specific 
events,  through  a  correlation  of  the  event  identification  given  in  the  Helsinki  Bulletin  (i.e., 
association  with  a  particular  mine)  and  that  resulting  from  a  cluster  analysis  (Riviere-Bar- 
bier,  1993).  A  result  of  this  study  is  a  set  of  well-controlled  events  that  will  further  con¬ 
strain  estimates  of  azimuth  bias. 

Mapping  observations  of  azimuth  bias  and  scatter  as  demonstrated  here  makes  it 
possible  to  build  azimuth-correction  tables  and  weight  tables  that  can  be  used  in  the  IMS 
to  improve  automated  locations. 
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Optimization  of  Polarization  Analysis  at  the  IMS/IRIS 
Stations  in  the  Former  Soviet  Union 

Anne  Suteau-Henson 


Introduction 

This  study  is  part  of  an  effort  to  integrate  data  from  six  three -component  (3-C) 
IRIS/IDA  stations  in  the  former  Soviet  Union,  AAK,  ARU,  GAR,  KTV,  OBN  and  TLY, 
within  the  automated  processing  of  the  Intelligent  Monitoring  System  (IMS,  Bache  et  al„ 
1990).  This  contribution  consisted  of  optimizing  the  measurement  of  polarization 
attributes  at  each  station,  for  use  in  phase  identification  and  azimuth  estimation.  In  a  previ¬ 
ous  study  (Riviere-Barbier  et  al.,  1992),  the  capabilities  for  phase  identification  were 
investigated  at  four  of  the  stations  (ARU,  GAR,  KTV,  and  OBN),  using  polarization  mea¬ 
surements  obtained  with  one  “recipe”  (i.e.,  set  of  processing  parameters)  for  all  stations. 
This  recipe  was  adapted  from  that  used  in  the  IMS  for  the  Scandinavian  arrays.  Station- 
specific  recipes,  significantly  different  from  the  “standard”  recipe,  were  found  to  optimize 
azimuth  estimation. 

In  this  study,  optimization  of  phase  identification  at  each  station  is  addressed.  The 
parameter  space  for  polarization  analysis  is  explored,  and  the  performance  of  the  resulting 
recipes  is  compared.  A  similar  experiment  is  carried  out  for  azimuth  estimation.  Finally, 
the  possibility  of  a  combined  optimization  of  phase  identification  and  azimuth  estimation 
is  considered.  For  each  station,  one  single  recipe  can  be  selected  that  jointly  optimizes 
both,  thus  considerably  simplifying  the  automated  processing. 

The  data  set  consists  of  the  detections  from  automated  processing  of  continuous 
data  covering  the  period  July  6  -  25,  1991.  At  each  station,  the  detections  were  divided 
into  four  groups:  Regional  S  (“S”)»  Regional  P  (“F’),  Teleseismic  (“T”),  and  Noise  (“N”), 
based  on  the  identification  given  by  the  analyst  who  reviewed  the  data.  The  data  set  is 
described  in  Table  2.  Not  enough  data  were  available  at  OBN  for  this  study.  'At  ARU,  a 
large  noise  detection  group  was  reduced  to  a  size  of  60  to  make  it  commensurate  with  the 
other  groups  at  that  station  (the  statistical  method  used  to  classify  the  detections  works 
best  for  groups  of  similar  sizes).  Except  for  AAK  and  GAR,  the  statistical  results  must  be 
considered  as  preliminary  due  to  small  populations. 
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Table  2:  Number  of  Detections  in  Data  Set 


Station 

Regional  S 

Regional  P 

Teleseismic 

Noise 

All 

AAK 

261 

254 

156 

56 

727 

ARU 

35 

20 

123 

60 

238 

GAR 

364 

341 

149 

203 

1057 

KIV 

101 

79 

62 

47 

289 

TLY 

67 

52 

104 

54 

277 

Exploration  of  the  Parameter  Space 

For  each  detection  in  the  data  set,  polarization  analysis  was  performed  to  obtain 
polarization  measurements  for  a  set  of  recipes,  including  the  standard  recipe  used  for  the 
IMS  arrays.  The  following  processing  parameters  were  varied: 

•  the  frequency  bands  to  include  in  the  analysis; 

•  the  length  of  overlapping  time  windows  in  which  to  measure  polarization  near 
the  detection; 

•  the  amount  of  overlap  between  time  windows; 

•  the  length  of  the  waveform  segment  analyzed  around  the  detection; 

•  the  time  delay  between  the  beginning  of  the  segment  and  the  onset  of  the  detec¬ 
tion. 

Other  parameters  (length  and  positioning  of  pre-detection  noise  segment)  were 
shown  to  have  little  effect  in  previous  studies  (Riviere-Barbier  et  al.,  1992).  For  the  new 
recipes,  the  length  of  the  filter  taper  was  increased,  and  the  SNR  threshold  (determining 
the  frequency  bands  selected  for  the  analysis)  was  set  higher,  compared  to  their  values  in 
the  standard  recipe.  The  method  for  automated  polarization  analysis  is  described  in  greater 
detail  in  Suteau-Henson  (1991). 

In  order  to  compare  the  performance  of  various  recipes  for  phase  identification, 
classical  discriminant  analysis  was  performed  as  a  method  of  classification.  Polarization 
measurements  were  used  as  predictors,  along  with  dominant  frequency  and  signal-to-noise 
ratio  (SNR).  Two  problems  were  considered: 

•  classification  of  the  three  signal  groups  (S  /  P  /  T),  not  including  coda  detections; 

•  classification  of  all  detection  groups  (i.e.,  including  coda  detections  and  the  noise 
detection  group). 
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As  the  first  problem  is  less  computationally-in tensive,  it  could  be  investigated  with 
a  laiger  set  of  recipes.  Eleven  different  recipes  were  tried  and  are  described  in  Table  3. 


Table  3:  Description  of  Recipes 


Recipe 

Freq. 

Band 

(Hz) 

Window 

Length 

(sec) 

Overlap 

Segment 

Length 

(sec) 

Onset 

Delay 

(sec) 

Filter 

Taper 

(sec) 

SNR 

Thresh. 

S 

1-4 

2 

0.5 

8 

4 

0.5 

1.5 

Al 

1.5 

0.5 

5 

0.5 

5 

2 

A2 

0.5-4 

1.5 

0.5 

5 

0.5 

5 

2 

A3 

0.5-8 

1.5 

0.5 

5 

0.5 

5 

2 

Gl 

1-2 

1 

0.5 

3 

1 

5 

2 

G2 

1-4 

1 

0.5 

3 

1 

5 

2 

G3 

1-8 

1 

0.5 

3 

1 

5 

2 

K1 

0.5-2 

1.5 

0.75 

3 

0.5 

5 

2 

K2 

0.5-4 

1.5 

0.75 

3 

0.5 

5 

2 

K3 

0.5-8 

1.5 

0.75 

3 

0.5 

5 

2 

A4 

0.5-8 

1.5 

0.5 

5 

1 

5 

2 

Exploring  the  effect  of  frequency  band  was  emphasized.  A  previous  study  (Riviere-Bar- 
bier  et  al 1992)  showed  that  the  best  recipes  for  azimuth  estimation  did  not  include  any 
frequencies  above  2  Hz.  On  the  other  hand,  wide-band  recipes  perform  better  for  phase 
identification,  as  expected  since  higher  frequencies  afford  better  discrimination  between 
regional  and  teleseismic  P,  in  particular.  The  eleven  recipes  can  be  briefly  described  as  fol¬ 
lows: 

(8)  Standard  Recipe  (S), 

(9)  ARU  Azimuth  Recipe  (Al),  i.e.,  as  optimized  by  Riviere-Barbier  et  al. 

(1992)  for  station  ARU, 

(10)  Same  as  (2),  but  including  the  2-4  Hz  band  (A2), 

(11)  Same  as  (3),  but  including  the  4-8  Hz  band  (A3), 

(12)  GAR  Azimuth  Recipe  (Gl), 

(13)  Same  as  (5),  but  including  the  2-4  Hz  band  (G2), 

(14)  Same  as  (6),  but  including  the  4-8  Hz  band  (G3), 
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(15)  KIV  Azimuth  Recipe  (Kl), 

(16)  Same  as  (8),  but  including  the  2-4  Hz  band  (K2), 

(17)  Same  as  (9),  but  including  the  4-8  Hz  band  (K3), 

(18)  Same  as  (4),  but  with  different  positioning  of  the  waveform  segment  (A4). 

Note  that,  although  the  higher  frequency  band  included  is  4-8  Hz,  in  practice,  an 
anti-aliasing  filter  applied  to  the  data  strongly  attenuates  frequencies  above  5  Hz. 

Performance  Comparison  for  Phase  Identification 

S/  PIT  Classification 

The  first  experiment  attempted  to  separate  regional  S,  regional  P,  and  teleseismic 
detections,  using  polarization  measurements  and  dominant  frequency  as  predictors.  The 
results  are  shown  in  Figure  7,  in  terms  of  overall  performance  at  each  station  for  each  rec¬ 
ipe.  The  following  observations  can  be  made: 

•  The  performance  for  phase  identification  of  the  recipes  that  were  previously 
found  to  optimize  azimuth  estimation  for  ARU,  GAR,  and  KIV  (“azimuth  reci¬ 
pes”,  from  Riviere-Barbier  et  al„  1992)  is  poor,  especially  at  GAR  and  KTV,  con¬ 
firming  our  previous  results. 

•  At  all  stations  except  ARU,  where  the  success  rate  is  low  regardless  of  recipe, 
adding  the  2-4  Hz  band  to  the  analysis  improves  the  performance. 

•  At  KTV  and  TLY,  results  are  further  improved  by  adding  the  4-8  Hz  band. 

•  The  “GAR  recipes”  (Gl,  G2,  G3)  that  do  not  include  the  0.5-1  Hz  “teleseismic” 
band,  perform  poorly  at  AAK  and  GAR. 

•  The  performance  of  the  standard  recipe  S  (that  includes  frequencies  from  1  to  4 
Hz)  is  generally  average,  but  poor  at  TLY. 

The  resulting  preferred  recipes  for  each  station  are  all  wide-band  and  include  fre¬ 
quencies  up  to  8  Hz.  At  a  given  station,  several  recipes  often  had  similar  overall  perfor¬ 
mance.  In  such  a  case,  their  relative  performance  for  the  different  phase  groups  was 
considered.  The  recipe  that  performed  best  for  regional  P  was  preferred.  For  example,  Fig¬ 
ure  8  shows  the  performance  at  ARU  for  each  phase  group  and  each  recipe,  and  illustrates 
the  selection  of  recipe  A3  rather  than  K3  at  this  station. 

S  /  P  /TIN  Classification 

Coda  detections  were  included  in  each  of  the  three  phase  groups,  and  the  noise 
detection  group  was  added  to  the  data  set  in  a  more  realistic  test  of  the  polarization  param¬ 
eters  for  phase  identification.  Also,  SNR  was  included  in  the  set  of  predictors.  At  each  sta- 
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Figure  7:  The  phase  identification  performance  (percentage  of  correctly  identified 
detections)  is  displayed  for  eleven  polarization  recipes  (described  in  the  text  and 
Table  3)  at  each  of  five  IMS/IRIS  stations.  Detections  in  the  three  groups  regional 
S,  regional  P,  and  tefeseismic  were  classified.  For  each  station,  the  standard  rec¬ 
ipe  (S),  the  optimum  recipe  for  azimuth  estimation  from  Riviere-Barbier  et  at. 
(1992),  when  it  applies,  and  the  preferred  recipe  from  this  experiment  are  high¬ 
lighted. 
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Figure  8:  The  phase  identification  performance  (percentage  of  correctly  identified  detections)  is  dsplayed  for 
eleven  polarization  recipes  (described  in  the  text  and  Table  3)  at  station  ARU.  Detections  in  the  three  phase 
groups,  regional  S,  regional  P,  and  Teleseismic,  were  classified.  In  spite  of  somewhat  lower  overall  performance 
(see  Figure  7),  recipe  A3  is  preferred  to  K3  due  to  its  higher  success  rate  for  the  regional  P  group. 


tion,  the  performance  for  each  of  the  four  detection  groups,  as  well  as  the  overall 
performance,  were  compared  for  three  recipes:  the  standard  recipe  (S),  the  azimuth  recipe 
(A)  when  available  for  this  station,  and  the  “best”  recipe  (B)  obtained  in  the  previous 
experiment  The  results  are  displayed  in  Figure  9  and  show  the  following: 

•  The  low-frequency  azimuth  recipes  (A)  have  a  low  overall  classification  rate  as 
in  the  previous  test  along  with  a  poor  performance  for  regional  P. 

•  The  high-frequency  recipe  B  outperforms  the  standard  recipe  S  at  TLY,  and  has  a 
comparable  overall  performance  at  the  other  stations. 

•  At  all  stations  except  ARU,  the  B  recipes  classify  regional  P  detections  better 
than  the  standard  recipe. 

Therefore,  at  each  station  except  ARU,  the  “best”  recipe  from  the  previous  exper¬ 
iment  still  qualifies  as  the  preferred  recipe  for  phase  identification  when  all  detections  are 
included  in  the  test  Its  advantage  over  the  standard  recipe  is  particularly  evident  at  TLY. 

Performance  Comparison  for  Azimuth  Estimation 

The  performance  for  azimuth  estimation  of  the  set  of  eleven  recipes  was  compared 
at  each  station.  The  data  sets  used  for  this  study  were  obtained  by  matching  the  event  ori¬ 
gins  after  analyst  review  with  origins  from  the  NEIC  Bulletin.  Azimuth  residuals 
(observed  azimuth  from  polarization  minus  azimuth  from  the  NEIC  location)  were  calcu¬ 
lated  for  teleseismic  P  arrivals.  Figure  10  describes  the  distribution  of  azimuth  residuals  at 
AAK,  ARU,  GAR,  and  TLY  for  each  recipe.  At  KIV,  very  poor  results  were  obtained 
regardless  of  recipe,  because  of  the  assumption  of  rectilinear  (instead  of  elliptical)  particle 
motion  used  in  the  polarization  analysis  (see  Riviere-Barbier  et  al.,  1992).  Therefore,  sta¬ 
tion  KIV  was  excluded  from  this  experiment.  The  figure  shows  the  following: 

•  The  standard  recipe  S  and  the  GAR  recipes  (Gl,  G2,  G3)  perform  relatively 
poorly  for  azimuth  estimation,  especially  at  AAK  and  TLY.  This  indicates  a  sen¬ 
sitivity  to  the  length  of  overlapping  time  windows  used  to  make  polarization 
measurements  (see  Table  3). 

•  The  azimuth  recipes  (Al  at  ARU  and  Gl  at  GAR)  are  not  optimum,  especially 
for  GAR. 

•  At  all  stations,  the  minimum  variance  is  obtained  with  a  recipe  including  the  2-4 
Hz  band,  and  at  AAK  and  GAR  the  best  recipe  includes  the  4-8  Hz  band. 

•  Except  for  the  poorly  performing  GAR  recipes,  the  variance  does  not  signifi¬ 
cantly  increase  as  higher  frequencies  are  included. 

Therefore,  based  on  these  results,  an  optimization  scheme  for  both  phase  identifi¬ 
cation  and  azimuth  estimation  utilizing  only  one  polarization  recipe  appears  feasible. 
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Figure  9:  The  phase  identification  performance  (percentage  of  correctly  identified  detections)  is  displayed  for 
each  of  the  four  groups  (regional  S,  regional  P,  Teleseismic,  and  Noise),  and  their  average,  at  each  of  five  IMS/ 
IRIS  stations.  The  results  are  shown  for  the  standard  recipe  “S’  the  optimum  recipe  “A"  for  azimuth  estimation 
(from  Riviere-Barbier  &  al.,  1992)  when  it  applies,  and  the  best  recipe  mBm  obtained  fortheS/P/T  classification 
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Figure  10:  The  distributions  of  azimuth  residuals  (azimuth  from  polarization  minus  azimuth  from  NEIC  location)  of 
teleseismic  P  arrivals  are  displayed  for  eleven  polarization  recipes  (described  in  the  text  and  Table  3),  at  each  of 
four  IMS/IRIS  stations. 


Combined  Optimization 

A  weighted  average  was  used  to  estimate  the  overall  performance  of  various  reci¬ 
pes  for  phase  identification  and  azimuth  estimation  combined.  At  the  3-C  IMS/IRIS  sta¬ 
tions,  phase  identification  is  more  important  than  azimuth  estimation,  and  therefore,  was 
given  a  weight  of  2.  Also,  as  the  correct  identification  of  regional  P  is  more  crucial  than 
that  of  other  detections,  it  was  assigned  a  weight  of  1.75. 

The  statistics  for  azimuth  estimation  displayed  in  Figure  10  were  converted  to 
measures  of  performance  (in  percent),  comparable  to  the  classification  performance  for 
phase  identification.  This  was  achieved  by  subtracting  the  lower  from  the  upper  95%  con¬ 
fidence  limits,  and  linearly  transforming  the  result  (x)  to  obtain  an  “azimuth  estimation 
performance”  y  (using  the  formula:  y  =  100  - 1.4  x).  The  results  of  combining  the  perfor¬ 
mance  for  phase  identification  with  that  for  azimuth  estimation  are  displayed  in  Figure  11 
for  the  set  of  eleven  recipes,  and  for  the  S  /  P  /  T  classification  experiment  described  in  the 
previous  section.  They  can  be  summarized  as  follows: 

•  The  GAR  recipes  (Gl,  G2,  G3)  and  the  standard  recipe  S  (except  at  ARU)  have  a 
poor  combined  performance. 

•  ARU  or  KIV  recipes  with  frequency  ranging  at  least  up  to  4  Hz  (i.e,  A2,  A3,  K2, 

K3)  perform  best 

Finally,  for  each  station,  a  subset  of  recipes  was  selected  and  tried  on  the  data  set 
including  the  coda  and  noise  detections  which  is  more  realistic  for  phase  identification  in 
the  IMS.  This  subset  included  the  standard  recipe  S,  the  azimuth  recipe  (at  stations  where 
it  applied),  the  best  recipe  from  phase  identification  optimization  alone,  and  any  other 
promising  recipe  based  on  the  results  of  the  earlier  test  of  combined  optimization  on  a 
smaller  data  set  (see  Figure  11).  Figure  12  shows  the  resulting  combined  performance.  In 
summary, 

•  at  AAK,  GAR  and  TLY  the  best  performance  is  obtained  for  wide-band  recipes, 
including  frequencies  from  0.5  to  8  Hz; 

•  at  ARU  a  better  performance  results  from  the  use  of  a  recipe  that  does  not 
include  the  higher  frequencies  (the  4-8  Hz  band). 

Therefore,  only  at  ARU  does  including  azimuth  estimation  in  the  optimization  sig¬ 
nificantly  alter  the  results.  This  may  be  due  to  the  small  number  of  regional  arrivals  in  the 
data  set  for  that  station,  and  a  larger  data  set  should  be  analyzed  to  reach  definite  conclu¬ 
sions. 
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Figure  11:  The  combined  performance  for  phase  identification  and  azimuth  esti¬ 
mation  is  displayed  for  a  set  of  eleven  polarization  recipes  (described  in  the  text 
and  Table  3),  at  four  IMS/IRIS  stations.  The  classification  problem  considered  is 
the  separation  of  regional  S,  regional  P,  and  teleseismic  detections. 
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Figure  12:  The  combined  performance  for  phase  identification  and  azimuth  esti¬ 
mation  is  displayed  for  four  of  the  IMS/IRIS  stations.  The  data  set  used  for  phase 
identification  included  all  four  detection  groups  (regional  S,  regional  P,  Teleseis- 
mic,  and  Noise).  At  each  station,  the  performance  was  estimated  for  the  standard 
recipe  S,  the  optimum  recipe  for  azimuth  estimation  from  Riviere-Barbier  et  al. 
(1992),  when  it  applied,  and  a  selected  subset  of  recipes  that  performed  well  for 
the  S/P /T  combined  optimization  experiment  (see  Figure  11).  The  resulting  pre¬ 
ferred  recipes  are  highlighted. 
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Conclusions 

Polarization  analysis  was  optimized  at  IMS/IRIS  3-C  stations  to  improve  its  use¬ 
fulness  in  the  IMS  automated  processing.  The  first  test  consisted  of  optimizing  the  recipe 
for  initial  phase  identification  alone.  The  frequency  band  used  in  the  analysis  was  found  to 
be  the  main  processing  parameter  affecting  the  performance.  The  recipes  previously  found 
to  optimize  azimuth  estimation  (Riviere-Barbier  et  al.,  1992),  as  they  do  not  include  fre¬ 
quencies  above  2  Hz,  performed  poorly,  especially  for  the  important  identification  of 
regional  P  detections.  Including  die  2-4  Hz  band  generally  increased  the  performance.  At 
TLY,  adding  the  4-8  Hz  band  (not  part  of  the  standard  IMS  recipe)  resulted  in  further  sig¬ 
nificant  improvement 

Next  optimization  of  azimuth  estimation  was  attempted  at  AAK,  ARU,  GAR,  and 
TLY.  This  time,  temporal  parameters  in  the  polarization  recipe  (such  as  length  of  overlap¬ 
ping  time  windows)  affected  the  performance  more  than  frequency  band.  The  recipes  that 
performed  best  included  higher  frequencies  (at  least  the  2-4  Hz  band)  as  for  phase  identifi¬ 
cation.  These  observations  indicated  that  at  these  stations  at  least  including  higher  fre¬ 
quencies  in  the  recipe  (which  improves  phase  identification)  did  not  significantly  degrade 
the  performance  for  azimuth  estimation. 

At  AAK,  GAR  and  TLY,  combined  optimization  of  phase  identification  and  azi¬ 
muth  estimation  was  achieved  with  a  wide-band  recipe,  including  frequencies  from  0.5  to 
8  Hz,  while  at  ARU,  a  recipe  with  a  narrower  band  (up  to  4  Hz  only)  gave  better  results. 
The  standard  IMS  recipe  performed  poorly,  mainly  due  to  its  inadequacy  for  azimuth  esti¬ 
mation,  which  may  be  related  to  its  lack  of  time  resolution.  No  single  recipe  was  found  to 
be  optimum  for  all  stations  considered.  This  is  not  surprising  as  the  signal  and  noise  char¬ 
acteristics  are  site-dependent  The  fact  that  at  a  given  station,  one  recipe  could  be  found 
satisfactory  for  both  phase  identification  and  azimuth  estimation  is  likely  due  to  the  fact 
that  they  both  rely  on  optimum  measurements  of  polarization.  The  polarization  algorithm 
insures  that  only  frequency  bands  with  SNR  above  a  given  threshold  contribute  to  the 
measurements,  so  that  higher  frequencies  do  not  degrade  azimuth  estimation  for  teleseis- 
mic  P,  for  example,  as  long  as  the  SNR  threshold  is  set  high  enough.  Based  on  these 
results,  it  appears  that  careful  parameterization  of  polarization  analysis  will  make  such  a 
combined  optimization  feasible  at  other  3-C  stations  as  well. 
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Stacking  of  Waveforms  for  Depth  Estimation 

Hans  Israelsson 

Introduction 

In  this  note  we  define  a  procedure  to  stack  waveforms  recorded  by  a  network  of 
seismic  stations  for  detection  of  depth  phases,  sP  and  pP.  Procedures  for  stacking  wave¬ 
forms  to  enhance  surface  reflected  phases  have  been  developed  and  tested  by  Roy  (1984), 
who  identified  depth  phases  from  the  maxima  of  stacked  signal  envelopes.  In  this  note  we 
aim  at  a  procedure  that  automatically  searches  for  a  depth  on  the  basis  of  surface  reflected 
phases.  Such  a  depth,  hereafter  called  a  ‘candidate  depth’,  is  given  as  an  interval  over 
which  possible  depth  phases  are  enhanced.  The  procedure  only  requires  die  epicenter  of 
the  event  to  be  processed,  and  the  first  arrival  times  of  the  associated  waveforms.  It  uses 
waveforms  recorded  in  the  teleseismic  window,  ie.,  30-90  degrees,  as  records  obtained  at 
shorter  distances  may  be  corrupted  by  die  complexity  of  wave  propagation  in  the  regional 
distance  range.  The  stacking  procedure  is  applied  to  data  from  die  world  wide  seismic  sta¬ 
tion  network  of  the  GSETT-2.  During  this  experiment,  as  described  by  Bratt  (1992),  wave¬ 
forms  recorded  at  a  world  wide  network  of  seismological  stations  were  collected  and 
compiled  from  April  22  through  June  2, 1991. 

Data 


For  the  analyses  here  we  selected  waveforms  that  could  be  associated  with  events 
reported  by  NEIS  with  a  body  wave  magnitude,  mb(NEIS)  >  4.0,  and  that  were  based  on 
observations  at  more  than  10  stations.  Furthermore,  events  with  depths  estimated  by  NEIS 
to  be  deeper  than  200  km  were  omitted. 

The  waveform  data  were  automatically  retrieved  from  the  GSETT-2  database  at 
the  Center  for  Seismic  Studies.  One  or  more  waveforms  in  the  teleseismic  window  were 
available  for  488  events,  the  geographical  distribution  of  which  is  shown  in  the  map  of 
Figure  13.  Two  of  the  events  were  underground  nuclear  explosions,  both  with  locations  in 
the  South  Pacific  and  marked  with  a  cross  in  Figure  13.  The  magnitude  of  the  events  were 
between  mb=4.1  and  6.4  and  had  depths  estimated  or  constrained  by  NEIS  between  5  and 
196  km,  apart  from  the  two  explosions.  In  all,  3,512  waveforms  were  processed.  They 
were  recorded  at  40  stations,  30  single  stations  and  10  arrays  with  apertures  varying  from 
a  few  to  a  few  tenths  of  km.  The  geographical  distribution  of  the  stations,  which  are  indi¬ 
cated  by  the  designated  station  codes  in  Figure  13,  is  highly  non-uniform.  In  fact,  almost 
50%  of  the  stations  are  concentrated  to  Europe  over  an  area  less  than  about  10%  of  the 
earth’s  surface.  The  distributions  of  waveforms  over  events  and  over  stations  were  also 
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GSETT-2  :  EVENTS  AND  STATIONS 
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Figure  13:  The  map  shows  the  geographical  distribution  of  the  events  and  the 
seismic  stations  for  which  waveform  data  were  analyzed.  The  two  under¬ 
ground  nudear  explosions  in  the  South  Padfic  are  marked  with  crosses. 


uneven.  The  number  of  waveforms  per  event  varied  between  1  and  34  with  a  median  of  6, 
and  the  number  of  waveforms  per  station  varied  between  2  and  280  with  a  median  of  78. 
Only  vertical  component  recordings  were  used  in  the  processing  and  analysis,  and  for 
array  stations  phased  beams  woe  used  to  the  largest  extent  possible. 

Depth  Traces 

Before  stacking  waveforms  associated  with  a  given  event,  each  waveform  was 
processed  in  several  steps  in  order  to  increase  the  signal-to- noise  ratio,  compensate  for  dif¬ 
ferences  in  delay  or  move  out  of  the  depth  phases,  and  to  minimize  the  possibility  of  mis¬ 
taking  a  PcP  phase  for  a  depth  phase.  Furthermore,  each  waveform  was  transformed  into  a 
binary  trace  with  a  STA/LTA  detector  to  enhance  possible  secondary  depth  phases  and  to 
avoid  destructive  interference  in  the  stacking  procedure  from  large  coda  levels  at  stations 
with  no  depth  phases.  The  steps  of  the  signal  processing  are  illustrated  in  Figure  14. 

Frequency  filtering 

It  has  been  observed  that  depth  phases  from  earthquakes  at  depths  below  200  km 
appear  to  have  less  high  frequency  than  the  direct  P  because  their  paths  traverse  the  struc¬ 
ture  above  the  focus  of  an  earthquake  twice  and  are  subject  to  more  significant  attenuation. 
(Blandford,  private  communication)  To  examine  similar  effects  for  earthquakes  at  shal¬ 
lower  depths.  Figure  15  shows  the  ratio  of  the  periods  of  pP/P  and  sP/P  for  phases  that 
were  positively  identified  during  the  GSETT-2.  The  ratios  are  plotted  as  a  function  of  the 
delay  time,  Lc.  arrival  times  of  pP  -  P  and  sP-P  respectively.  The  ratios  have  a  large  scatter 
and  there  is  no  consistent  difference  in  the  periods  of  pP  and  sP  on  one  hand  and  of  direct 
P  on  the  other.  However,  the  smoothed  trends  obtained  by  robust  regression  (Cleveland, 
1979)  suggest  that  the  pP  and  sP  appear  to  have  slightly  more  low  frequency  than  does 
direct  P. 

In  order  to  maintain  as  broad  a  frequency  range  as  possible  and  to  allow  for  differ¬ 
ences  in  noise  and  receiver  responses  among  the  stations,  the  limits  of  the  band  pass  filter 
(3  pole  phase  free  Butterworth)  were  based  on  spectral  analysis  of  noise  prior  to  the  onset 
of  the  initial  P. 

Binary  trace 

Stacking  aimed  at  enhancing  depth  phases  that  are  based  on  bandpass  filtered 
waveforms  or  waveform  envelopes  may  suffer  from  destructive  interference  from  wave¬ 
forms  with  high  coda  levels  and  with  no  clearly  developed  surface  reflected  phase.  There¬ 
fore  we  applied  a  STA/LTA  detector,  with  windows  of  1  and  5  s  length  respectively,  to  the 
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Figure  14:  The  signal  processing  and  transformation  in  the  time  domain, 
upper  frame,  and  in  the  depth  domain,  lower  frame.  The  time  domain  process¬ 
ing  involves  bandpass  filtering  to  enhance  the  signal-to-noise  ratio,  forming  of 
a  STAA.TA  trace,  masking  of  window  for  expected  arrival  time  of  PcP,  and 
transformation  into  a  binary  trace.  The  binary  trace  in  the  time  domain  is 
mapped  into  two  binary  traces  in  the  depth  domain,  one  each  for  pP  and  sP. 
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Figure  15:  The  scatter  diagrams  shows  the  ratios  of  the  periods  ofpP/P, 
to  the  left,  and  of  sP/P  to  the  right,  as  a  function  of  the  delay  times. 
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waveforms.  This  STA/LTA  trace  is  then  converted  into  a  binary  trace,  which  is  0  or  1  if  the 
STA/LTA  ratio  is  less  than  or  greater  than  a  threshold  of  1.75,  respectively.  In  Figure  16 
we  have  compiled  all  the  STA/LTA  traces  into  a  record  section  as  a  function  of  epicentral 
distance.  The  part  of  the  traces  with  STA/LTA  greater  than  the  threshold  are  indicated  as 
black  dots.  The  expected  arrival  time  of  PcP  for  depths  between  0  and  200  km  is  indicated. 
No  clear  clustering  of  detected  phases  appears  along  the  lines  indicated  for  the  expected 
PcP,  and  the  section  gives  a  general  impression  of  randomness  with  no  other  significant 
features.  The  generally  white  area  5-10  s  following  the  direct  P,  may  be  partly  due  to  the 
saturation  of  the  STA/LTA  ratio  caused  by  the  primary  phase  which  usually  has  the  largest 
amplitudes. 

Masking  of  PcP  window 

Even  if  the  PcP  phase  does  not  stand  out  as  a  striking  feature  on  the  record  section 
in  Figure  16,  we  mask,  or  zero  the  binary  trace  for  a  2  s  window  centered  at  the  expected 
arrival  time  of  PcP.  This  avoids  mis-identifying  PcP  phases  as  depth  phases  in  cases  where 
PcP  is  clearly  developed. 

Transformation  into  depth  traces 

Before  binary  traces  calculated  from  waveforms  recorded  at  different  distance 
ranges  could  be  stacked  together  we  made  a  transformation  from  the  time  domain  to  the 
depth  domain.  This  compensated  for  possible  differences  in  delay  or  move-out  of  the 
depth  phases  relative  to  the  direct  P  for  different  epicentral  distances.  The  move-out  also 
strongly  depends  on  the  depth  of  the  event  For  example,  pP  move-out  is  about  1  and  9  s 
between  30  and  90  degrees  for  earthquakes  at  50  and  200  km  depth,  respectively. 

Each  binary  trace  was  thus  mapped  into  two  depth  traces  one  for  the  pP  and  one 
for  sP.  A  depth  trace  for  pP,  for  example,  is  obtained  by  converting  the  delay  time  from  the 
onset  of  direct  P  into  the  corresponding  depth  using  a  travel  time  table.  For  this  study  we 
used  the  JB  travel  time  tables.  As  the  synthetic  example  in  Figure  14  illustrates  for  a  case 
with  clearly  developed  pP  and  sP  phases,  the  binary  depth  traces  for  pP  and  sP  will  have 
signals  that  line  up  at  the  same  depth,  around  70  km  in  this  example. 

Stacking  of  Binary  Traces  and  Candidate  Depths 

After  calculating  the  binary  traces,  they  were  combined  into  stacks  for  pP  and  sP 
by  straight  summation.  Only  binary  traces  corresponding  to  waveforms  with  signal-to- 
noise  ratios  larger  than  3.0  were  used  in  the  stacking.  Selected  examples  of  binary  and 
stacked  traces  together  with  the  recorded  waveforms  are  shown  in  Figure  17  a-c  for  three 
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Figure  16:  STA/LTA  traces  plotted  as  a  function  of  epicentral  distance  for  all  the  waveforms.  Only 
the  parts  of  the  traces  with  the  STA/LTA  above  the  threshold,  1.75,  were  marked  as  black  dots. 


JSSyNY^vf-*  *ft*M lW«W****U*..tf 


mm  mnu 


Figure  17:  The  figure  shows  waveforms  and  binary  traces  sections  for  three  events  in  the 
Afghanistan-USSR  border  region  on  :  (a)  91/05/13  at  23:55:33.7,  (b)  91/04/26  at 
22:24:05.3,  and  (c)  91/05/13  08:32:43.8.  Station  codes  and  epicentral  distances  are  indi¬ 
cated  to  the  right  of  each  waveform.  The  waveform  sections  also  include,  at  the  top,  enve¬ 
lope  stacks  of  pP  and  sP  as  a  function  of  depth  and  the  standard  deviation  of  the  first 
arrival  time  P  residuals  as  a  function  of  depth  (SIGMA  P  RES  in  the  diagrams).  The  signal- 
to-noise  ratio,  SNR,  and  the  P  residual  for  the  NEIS  hypocenter  solution  are  shown  to  the 
right  for  each  waveform. 
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earthquakes  near  the  Afghanistan-USSR  border  region.  The  three  figures  have  the  same 
format  and  each  one  includes  one  section  with  the  recorded  waveforms  (bandpass  filtered) 
and  one  section  with  the  binary  traces.  In  Figure  17a  for  the  event  on  May  13  at 
23:55:33.7,  clear  depth  phases  can  be  seen  on  many  of  the  filtered  traces.  The  signal-to- 
noise  ratio  and  arrival  time  residual(s)  are  indicated  at  the  end  of  each  trace.  As  wave¬ 
forms  with  SNR  less  than  3  were  not  used  in  the  stacking,  the  stacked  traces  are  based  on 
a  total  of  10  traces,  6  of  which  are  in  Fennoscandia,  all  with  quite  similar  appearance. 
Above  these  traces  are  also  stacks  of  the  envelopes  of  the  waveforms  transformed  to  a 
depth  scale  (also  indicated  at  the  top).  There  is  one  stack,  respectively,  for  pP  (thick  heavy 
line)  and  sP  (thick  shaded  line).  The  two  stacks  have  two  clear  secondary  peaks  that  line 
up  at  a  depth  of  about  25  km.  The  binary  trace  section  of  Figure  17a  also  shows  clear 
detections  of  the  depth  phases  on  the  individual  station  traces  as  well  as  on  the  binary 
stacks  at  the  top  of  the  section.  In  the  binary  trace  section,  the  horizontal  scale  for  all 
traces  is  depth.  Figure  17b  is  almost  the  opposite  of  Figure  17a  with  regard  to  observabil¬ 
ity  and  detectability  of  depth  phases.  Although  this  is  a  somewhat  larger  event  with  more 
waveforms  that  generally  also  have  larger  SNR,  any  clearly  identifiable  depth  phases  that 
apr«*ar  on  most  of  the  traces  can  not  be  identified.  13  of  the  17  traces  with  sufficient  SNR 
xcorded  at  European  stations  and  they  all  have  similar  features.  Possible  depth 
pii^es  can  only  be  seen  at  three  stations  (HFS,  GEC2  and  KSP).  The  stacked  envelope 
shows  a  maximum  at  about  10  km  depth  for  pP,  but  there  is  no  peak  on  the  sP  stack  that 
lines  up  with  this  maximum.  On  the  binary  trace  section  few  secondary  phases  were 
picked  by  the  STA/LTA  detector  and  although  they  line  up  to  show  a  secondary  maximum 
on  the  stack,  the  width  of  the  initial  phase  is  too  wide  to  clearly  separate  it.  Finally,  the 
example  in  Figure  17c  is  the  most  typical  of  the  data  analyzed  here  with  regard  to  the 
number  of  available  waveforms.  Secondary  phases  can  be  noted  on  at  least  two  of  the  five 
filtered  records  and  they  add  up  constructively  to  show  a  maximum  in  the  depth  envelope 
stack  corresponding  to  a  depth  of  about  100  km  from  pP.  However,  this  maximum  is 
somewhat  ambiguous  as  there  is  a  maximum  around  45  km  as  well.  The  stacking  on  the 
binary  section  separates  clearly  between  these  two  possibilities.  Only  the  depth  around 
100  km  is  now  apparent  Clearly  this  in  not  sufficient  to  positively  estimate  the  depth.  As 
the  example  in  Figure  17c  is  the  most  representative  of  a  typical  event  analyzed  in  this 
study,  we  use  the  stacked  binary  depth  trace  to  identify  ‘candidate’  depth  intervals  as 
described  in  the  following,  rather  than  attempting  to  assign  an  exact  depth  to  each  event 

Such  candidate  depths  were  defined  from  the  stacked  binary  traces  with  clearly 
defined  signals  (see  Figure  18).  For  the  purpose  of  the  analysis  here,  we  used  the  follow¬ 
ing  automatic  procedure.  A  given  stacked  depth  trace,  for  pP  or  sP,  was  searched  for  a 
unique  signal  maximum  after  the  signal  corresponding  to  the  initial  P  had  come  to  zero. 
That  is  to  say,  an  uninterrupted  depth  range  with  a  maximum  larger  than  any  other  value  of 


45 


the  depth  trace,  except  for  the  initial  P.  The  depths  corresponding  to  the  beginning  and  end 
of  the  interval  in  which  the  maximum  occurred  define  a  depth  interval  for  the  trace.  As  the 
pP  often  is  recorded  more  clearly  than  sP  phases,  the  candidate  depths  from  the  sP  trace  is 
just  a  non-linear  translation  of  the  depth  of  pP  trace. 

Applicability 

In  order  to  get  some  idea  of  the  performance  of  the  stacking  procedure  in  defining 
‘candidate’  depths,  we  examined  the  events  for  which  depths  could  and  could  not  be 
obtained  without  any  consideration  of  the  quality  of  the  obtained  candidate  depths  them¬ 
selves.  The  applicability  may  be  determined  by  several  factors,  such  as  the  actual  depth  of 
the  event,  the  number  of  available  waveforms,  and  the  distribution  of  the  stations  around 
the  source  azimuth.  Furthermore,  there  may  be  significant  regional  effects.  It  has,  for 
example,  been  suggested  that  scattering  and  defocusing  of  P  and  S  by  lateral  variations  in 
the  crust  and  upper  mantle  above  the  focus  may  cause  significant  reduction  in  the  ampli¬ 
tude  of  surface  reflections.  The  distribution  of  the  events  and  waveform  data  used  in  this 
study  was  not  sufficient  to  study  these  factors  in  isolation  or  to  apply  additive  models  that 
would  separate  the  effects  of  various  factors.  Therefore,  we  look  at  each  factor  separately. 

Figure  18a  shows  the  distribution  of  depth,  as  determined  by  NEIS,  for  the  events 
for  which  ‘candidate  depths’  could  and  could  not  be  defined.  Figure  18b  shows  the  per¬ 
centage  of  events  with  candidate  depths  as  a  function  of  depth.  The  data  in  these  two  fig¬ 
ures  are  limited  to  events  for  which  depth  was  determined  from  three  of  more  waveforms. 
We  note  that  most  of  the  NEIS  depths  were  constrained  to  either  10  or  33  km  (Figure  18a) 
and  that  a  significant  portion  of  the  events  that  are  constrained  had  ‘candidate’  depths.  The 
data  in  Figure  18b  does  not  suggest  any  clear  effect  of  focal  depth,  although  one  might 
have  expected  that  the  chances  of  defining  a  ‘candidate’  depth  would  increase  with 
increasing  depth. 

Figures  19a-b  show  the  geographical  distribution  of  the  events  with  and  without 
‘Candidate’  depth  marked  with  plus  signs  and  open  circles,  respectively.  For  continental 
events  in  Eurasia,  a  regional  pattern  is  evident.  For  example,  many  events  in  the  Hindu 
Kush  area  and  around  the  Indonesian  archipelago  have  candidate  depths,  whereas'few 
events  in  Southern  Iran  have  candidate  depths.  Candidate  depths  were  not  obtained  for  the 
two  underground  nuclear  explosions  in  the  South  Pacific. 

In  Figure  20  we  show  the  success  and  failure  rates  as  functions  of  the  number  of 
waveforms  and  azimuthal  coverage.  Azimuthal  coverage  is  defined  as  the  number  of  45 
degree  angular  sectors  around  the  source  to  station  azimuth  for  which  there  is  one  or  more 
stations.  Before  counting  the  number  of  sectors,  they  were  centered  at  the  median  azimuth 
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Figure  18:  (a)  The  figure  shows  how  the  candidate  depth  is  defined  from  a 
pP  depth  trace  stack. 
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Figure  18:  (b)  The  barplot  shows  the  number  of  events  with  and  without  can¬ 
didate  depth  as  a  function  of  depths  reported  by  NEIS.  Most  of  the  NEIS 
depths  were  constrained  to  33  km  and  a  significant  portion  of  those  had  candi¬ 
date  depths. 


48 


so 


too 


T 

lONomjoc 


Figure  19:  The  geographical  distribution  in  the  Western  (a)  and  Eastern  (b) 
hemisphseres  of  events  with  and  without  candidate  depths. 


of  the  stations.  This  means  that  cases  where  all  stations  are  clustered  within  45  degrees 
will  be  counted  as  one  45  degree  sector,  whereas  an  arbitrary  reference  of  the  azimuth  sec¬ 
tor  might  have  yielded  two  45  degree  sectors.  As  in  Figure  20  plus  signs  and  open  circles 
mark  cases  where  a  candidate  depth  could  and  could  not  be  obtained  respectively.  The 
success  rates,  defined  as  the  percent  of  events  for  which  ‘candidate’  depths  were  defined 
are  given  vertically  and  horizontally  as  a  function  of  the  number  of  45  degree  sectors  cov¬ 
ered  and  number  of  waveforms  available,  respectively.  With  four  or  more  sectors  the  suc¬ 
cess  rate  is  about  or  above  50%  and  9  or  more  waveforms  are  required  to  consistently 
achieve  this  rate.  We  note  from  the  plus  signs  and  open  circles  that  there  is  a  significant 
number  of  events  for  which  the  number  of  waveforms  is  smaller  than  9  but  the  number  of 
45  degree  sectors  covered  is  larger  than  4  for  which  candidate  depths  were  obtained. 

Comparison  of  Candidate  and  NEIS  Depths 

The  quality  of  the  derived  ‘candidate’  depths  are  difficult  to  assess  due  to  the  obvi¬ 
ous  lack  of  ground  truth  in  many  instances.  In  Figure  21  we  compare  the  NEIS  depths  and 
depths  taken  as  the  average  of  the  lower  and  upper  limits  of  the  candidate  depths.  Compar¬ 
isons  are  made  as  a  function  of  the  number  of  phases  in  the  stacked  trace,  1, 2, 3  and  more 
than  3.  The  number  of  events  for  each  case  is  indicated  in  the  diagrams.  Somewhat  arbi¬ 
trarily,  we  accepted  a  difference  of  20  km  between  the  two  types  of  depths  for  agreement 
(’’agree.”  in  the  diagrams).  This  large  acceptance  difference  was  partly  based  on  the  fact 
that  a  large  percentage  of  the  NEIS  depths  were  constrained  to  10  or  33  km.  The  compari¬ 
sons  in  Figure  21  show  reasonable  agreement  for  3  or  more  phases,  although  there  are  4 
events,  2  at  10  and  2  deeper  than  100  km  by  NEIS,  for  which  there  is  a  clear  discrepancy 
between  the  ‘candidate’  and  the  NEIS  depths.  The  agreement  rates  are,  as  one  would 
expect,  smaller  for  stacks  based  on  1  and  2  phases  only. 

Earthquake  Source  Radiation  Pattern  and  Depth  Phases 

As  has  been  noted  by  Douglas  et  al.  (1990),  theoretical  modelling  of  seismograms 
from  double  couple  sources  radiating  pulses  of  around  1  s  duration  show  that  the  most 
prominent  arrivals  on  short  period  seismograms  should  be  the  direct  P,  pP,  and  sP.  Indeed, 
Pearce  (1980)  developed  a  method  to  constrain  the  focal  mechanism  from  relative  ampli¬ 
tudes  of  P,  pP  and  sP.  However,  the  observations  in  the  past  of  the  character  of  P  seismo¬ 
grams  show  that  this  model  is  frequently  contradicted. 

Moment  tensors  have  been  determined  by  HVD  for  many  of  the  earthquakes  ana¬ 
lyzed  here  for  which  waveform  data  were  available  at  two  or  more  stations.  These  moment 
tensor  solutions,  which  are  obtained  from  long  period  data,  were  used  to  calculate  ‘theo¬ 
retical’  amplitude  ratios  for  pP/P  and  sP/P.  Algorithms  for  the  radiation  from  moment  ten- 
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Figure  20:  Success  rates,  given  as  percentages,  to  define  candidate  depths 
as  a  function  of  the  number  of  waveforms  and  azimuthal  coverage.  Azimuthal 
coverage  is  given  as  the  number  of  45  degree  sectors  around  the  source  to 
station  azimuth  for  which  there  is  one  or  more  station  with  waveform  data. 
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sors  by  Kennett  (1988)  together  with  the  procedure  to  calculate  the  reflection  coefficients 
of  surface  reflected  phases  defined  by  Pearce  (1980)  were  employed  for  this  purpose.  The 
theoretical  ratios  were  compared  with  those  obtained  from  the  waveform  data.  The 
'observed*  ratios  were  all  calculated  for  a  frequency  band  0.8-2.0  Hz.  In  the  calculations 
we  used  the  depths  reported  by  NEIS  and  omitted  events  with  constrained  depths,  i.e. 
depths  of  10  or  33  km. 

In  Figure  22  we  show  the  correlation  coefficient  of  calculated  and  observed  ratios 
plotted  as  a  function  of  die  difference  in  pP  delay  time  and  the  source  duration  as  reported. 
Source  duration  was  subtracted  from  the  pP  delay  because  there  may  be  interference 
between  direct  P  and  pP  for  events  at  very  shallow  focal  depth  and  with  long  source  dura¬ 
tions.  We  show  the  correlation  coefficient  of  pP  and  sP  ratios  separately  and  also  the  geo¬ 
metric  mean  of  these  two  correlation  coefficients.  Only  events  for  which  five  or  more 
waveforms  were  available  arc  included  in  the  comparison.  High  mean  correlation  coeffi¬ 
cients  were  observed  for  only  for  2  out  of  36  events.  Figure  23  shows  the  actual  observed 
and  calculated  ratios  for  one  of  the  events,  on  April  26  at  17:36:59.4, 7.4  N  127.7  E  and  74 
km  depth.  It  should  be  noted  that  14  of  the  events  occurred  at  such  a  depth  and  had  such 
duration  that  interference  between  P  and  pP  is  likely.  The  poor  correlation  between  the 
observed  and  calculated  ratios  could  be  due  to  several  factors,  but  the  results  illustrate  the 
unpredictability  of  depth  phases  from  simple  models  of  radiation  patterns. 

Concluding  Remarks 

In  this  note  we  have  defined  an  automatic  procedure  of  waveform  stacking  for  esti¬ 
mating  focal  depth  of  a  seismic  event.  As  input  it  requires  the  epicenter  and  a  preliminary 
estimate  of  the  focal  depth  together  with  recorded  waveforms  and  first  arrival  times.  The 
procedure  provides  a  'candidate  depth’,  defined  as  an  interval  over  which  stacked  binary 
traces  have  a  maximum  value.  The  procedure  should  be  seen  as  a  tool  assisting  analysts 
working  with  global  network  waveform  data  in  their  effort  to  determine  focal  depths  by 
presenting  possible  depth  solutions. 

The  stacking  procedure  was  tested  on  waveform  data  for  488  events  collected  dur¬ 
ing  the  GSETT-2.  The  capability  to  define  a  candidate  depth  depends  on  several  factors 
and  a  large  azimuthal  coverage  of  recording  stations  is  identified  as  one  such  factor  of 
importance.  A  candidate  depth  could  be  defined  for  about  50%  of  the  events  for  which  the 
azimuthal  coverage  was  180  degrees  or  more. 
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Figure  21:  Depths  taken  as  the  average  of  the  lower  and  upper  values  of  the 
candidate  depths  (‘pP  DEPTHS’  in  the  diagrams)  compared  with  NEIS  depths. 
Comparisons  are  made  as  a  function  of  the  number  of  phases  in  the  stacked 
trace,  1,  2,  3,  and  more  than  3.  Somewhat  arbitrarily,  a  difference  of  20  km 
between  the  two  types  of  depth  is  accepted  for  agreement  ('agree.' given  as 
percentages  in  the  diagrams). 
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Figure  22:  Comparison  of  theoretical  and  observed  amplitude  ratios  pP/P 
and  sP/P  for  the  event  on  91/04X6  1 7:36:59  with  latitude  7.40  N  and  iongr- 
tude  126. 71  E  and  mb(NEIS)=5.5. 
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Figure  23:  Correlation  coefficients  for  observed  and  theoretical  pP/P  and  sP/ 
P  amplitude  ratios  plotted  against  the  difference  between  pP  delay  time  and 
source  duration  for  GSETT-2  events  that  were  recorded  at  more  than  5  sta¬ 
tions. 


55 


The  capability  of  the  automatic  stacking  procedure  may  be  improved  by  tuning 
some  of  its  parameters  as  no  attempts  were  made  to  optimize  the  performance.  The  param¬ 
eter  settings  are  somewhat  arbitrary,  but  represent,  it  is  felt,  reasonable  values.  Further¬ 
more  the  procedure  can  be  extended  to  include  recordings  at  distances  closer  than  30 
degrees. 

Comparisons  of  observed  amplitude  ratios  for  pP/P  and  sP/P  phases  based  on  the 
amplitude  level  at  the  delay  time  predicted  from  focal  depths  determined  by  NEIS  with 
those  calculated  from  moment  tensor  solutions  determined  by  HVD  for  36  of  the  events 
showed  little  correlation.  Data  for  only  two  of  the  events  appeared  to  have  clear  correla¬ 
tions.  The  lack  of  correlation  could  be  attributed  to  several  factors;  incorrect  NEIS  depths; 
invalidity  of  the  extrapolation  of  the  moment  tensors  from  long  to  short  period  data;  and/ 
or  interference  of  direct  P  and  pP  due  to  shallow  focal  depth  and  long  source  duration. 
Whatever  the  reason  may  be,  the  results  illustrate  the  difficulty  of  reconciling  observed 
amplitude  ratios  from  those  calculated  from  the  radiation  pattern  of  the  source,  although 
for  some  events  there  is  apparent  agreement  between  observed  ratios  and  those  predicted 
from  simple  modelling. 
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Constraining  Origin  Time  with  S/P  Travel  Time  Ratios 
and  Depth  Estimation 

Hans  Israelsson 


Introduction 

In  this  note  we  study  with  synthetic  examples  the  constraint  on  focal  depth  provided  by 
origin  times  determined  from  travel  time  ratios  of  P  and  S  waves.  We  also  analyze  travel 
time  ratios  for  Pn  and  Sn  waves  recorded  by  the  IMS  arrays  and  determine  and  compare 
depths  for  about  70  events  -  presumed  earthquakes  as  well  as  presumed  chemical  explo¬ 
sions  in  and  near  Fennoscandia. 

Constraining  Origin  Time  With  S/P  Travel  Time  Ratio 

It  is  common  knowledge  that  accurate  depth  determinations  based  on  only  P  wave  arrival 
times  often  are  difficult  to  obtain  (Lomnitz,  1977).  This  is  due  to  the  limitations  of  such 
observations  to  separate  depth  and  origin  time,  which  sometimes  result  in  a  near  indeter¬ 
minacy  of  the  two  parameters.  This  can  be  illustrated  by  the  following  somewhat  simpli¬ 
fied  consideration.  Assume  that  the  epicenter  of  an  event,  the  depth  of  which  is  to  be 
determined,  is  known  and  that  n  observed  P  arrival  times,  to,  with  similar  standard  errors, 
s,  are  available  for  the  determination.  Forming  the  residuals,  r,-  =  to,-  -  to,-,  between  to,-  and 
theoretical  arrival  times,  to,-,  calculated  from  a  tentative  origin  time  and  depth  and  a  given 
travel  time  model,  the  adjustments  in  depth,  dh,  and  origin  time,  dt ;  are  related  to  the  resid¬ 
uals,  r,-,  by  a  standard  regression  relation:  r^dt  +  dh  (dT/dh)j.  The  depth  derivative  dTIdh 
is  obtained  from  the  travel  time  model  (derivative  for  the  ISAPEI91  model  is  shown 
graphically  in  Figure  24  as  a  function  of  distance  for  different  depths).  The  variance  of  the 
dh  estimated  from  this  regression  can  be  written  as: 


This  means  that  the  variance  of  the  depth  is  inversely  proportional  to  the  variance  of  the 
depth  travel  time  derivatives.  The  curves  in  Figure  1  show  that  the  derivatives  are  almost 
constant  with  little  variance  beyond  30  degrees  and  that  events  based  only  on  P  observa- 
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Figure  24:  The  derivative  dT/dh  as  a  function  of  distance  for  depths 
between  12  and  75  km. 
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tions  in  this  distance  range  clearly  will  have  depths  that  are  poorly  estimated. 

If  the  origin  time  can  be  determined  from  observations  of  phases  other  than  P,  however,  a 
location  constrained  to  that  origin  time  may  provide  an  accurate  estimate  of  the  focal 
depth  (Evemden,  1969).  For  example,  arrival  times  of  S  and  P  waves,  ts  and  tp,  can  be 
used  to  determine  the  origin  time,  T,  if  the  travel  time  ratio  for  S/P  waves,  r,  is  known:  T  = 
(r  ts  -  tpV(r-l).  Travel  time  ratios  for  S/P  waves  are  usually  quite  stable  over  some  distance 
range,  as  illustrated  by  the  calculated  ratios  for  the  travel  times  of  JB,  IASPEI91,  and  IMS 
in  Figure  25.  The  standard  deviations  for  the  ratios  vary  from  fractions  of  a  percent  to 
about  1%.  Clearly  lateral  variations  in  travel  time  will  increase  the  scatter,  but  travel  times 
can  be  calibrated  to  improve  precision. 

Use  of  S  arrival  times  by  standard  location  algorithms  will  automatically  improve  the  pos¬ 
sibility  of  separating  depth  and  origin  time.  Over  the  years,  several  rules  of  thumb  on  the 
sensitivity  of  focal  depth  solutions  to  the  inclusion  of  S  arrival  times  have  been  developed 
and  studied  (Gomberg  et  al.  1990).  From  the  non-linear  nature  of  the  hypocenter  determi¬ 
nation  problem,  however,  it  is  not  obvious  that  the  straightforward  and  standard  incorpora¬ 
tion  of  S  arrival  times  provides  depth  solutions  similar  or  equal  to  those  based  on  an  origin 
time  constrained  by  S/P  arrival  time  ratios. 

In  Figure  26  we  show  comparisons  based  on  synthetic  calculations,  which  suggest  that  a 
solution  based  on  a  constrained  origin  time  may  indeed  have  smaller  standard  error  than  a 
standard  solution,  at  least  under  certain  circumstances.  In  the  figure,  the  standard  devia¬ 
tions  of  estimated  depths  are  compared  for  two  hypothetical  networks,  both  with  5  stations 
at  about  the  same  distance  from  the  epicenter  (15  degrees).  The  depth  solutions  are  based 
on  P  arrival  times  at  all  stations  and  on  S  times  at  one  and  five  stations  in  the  upper  and 
lower  frames  respectively.  Algorithms  in  a  hypocenter  determination  program,  (LOCEX) 
originally  developed  by  Julian  (1974)  were  employed  in  the  simulations.  Gaussian  arrival 
time  errors  obtained  by  a  random  noise  generator  were  added  to  theoretical  arrival  times 
and  depths  were  estimated  for  3000  event  randomizations  at  depths  going  from  10  to  50 
km  in  steps  of  2  km.  In  the  simulations  we  assumed  that  the  S/P  travel  time  ratio  was  with¬ 
out  error.  The  dashed  lines  in  the  two  frames  to  the  right  represent  standard  error  equal  to 
actual  depth.  The  depths  obtained  from  constrained  origin  times  have  systematically 
smaller  standard  deviation,  up  to  about  a  factor  of  1 .5  -  2,  than  the  standard  depth  solu¬ 
tions.  The  differences  are  less  pronounced  in  the  case  with  5  Sn  phases.  Constraining  the 
origin  time  with  the  Sn/Pn  travel  time  ratio  in  the  case  with  only  one  Sn  phase  gives  much 
more  weight  to  the  Sn  observation  than  the  free  depth  location  does  and  it  may  therefore 
be  possible  to  reduce  the  standard  deviations  with  the  free  depth  solution  if  the  Sn  obser¬ 
vations  were  given  more  weight. 
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Figure  25:  The  ratio  for  travel  times  of  S  and  Pas  a  function  of  epicentral 
distances  for  different  velocity  models. 
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Figure  26:  Comparisons  of  the  standard  deviations  of  the  depth  estimates 
obtained  with  and  without  constrained  origin  times  for  two  synthetic  exam¬ 
ples,  with  recorded  phases  of  Pn  and  Sn  as  illustrated  by  the  maps  to  the 
left. 
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Ts/Tp  Ratios  From  IMS  Data 

In  order  to  get  some  idea  of  the  stability  of  the  S/P  travel  time  ratio  we  analyzed  arrival 
times  for  Pn  and  Sn  in  the  IMS  data  base.  The  results  are  summarized  in  Figure  27  by 
comparisons  of  origin  times  determined  from  Pn  and  Sn  arrival  times  at  the  three  arrays  in 
Fennoscandia.  For  each  pair  of  arrays,  a  ratio,  r,  was  fit  to  the  computed  origin  times  of 
common  events  that  minimized  the  differences  in  origin  times.  The  bias  and  the  standard 
deviations,  SD,  given  for  each  array  pair  in  the  figure  were  obtained  as  the  median  of  the 
differences  in  origin  time  and  a  robust  estimate  of  the  standard  error.  We  note  that  the  val¬ 
ues  vary  somewhat  among  array  pairs  with  a  range  of  about  1%.  Furthermore,  the  values 
that  involve  the  ARCESS  array  (ARAO)  are  also  slightly  higher  than  1.756,  which  is  the 
average  value  for  Sn/Pn  over  the  distance  range  5  to  20  degrees  in  the  travel  time  tables 
used  by  IMS.  The  estimates  of  the  S/P  travel  time  ratios  obtained  here  are  completely 
independent  of  the  IMS  or  any  other  travel  time  model.  The  estimated  standard  deviation 
of  the  origin  time  based  on  one  observation  is  2-2.5  s. 

A  large  standard  error  in  the  S/P  travel  time  ratio  will  clearly  result  in  a  large  standard 
error  of  the  origin  time  as  illustrated  by  Figure  28 ,  which  shows  curves  for  standard  error 
in  origin  time  as  a  function  of  epicentral  distance.  The  curves  were  obtained  by  simula¬ 
tions  using  a  Gaussian  random  number  generator  and  assuming  standard  errors  of  0.7  and 
1.4  s  in  P  and  S  arrival  times,  respectively. 

If  we  assume  that  the  standard  errors  in  the  arrival  times  of  P  and  S,  sp  and  ss  respectively, 
are  independent  and  that  ss  =  2  sp,  we  can  estimate  the  standard  error  in  sp  from  the  esti¬ 
mates  of  the  standard  error  of  the  origin  time,  sT,  in  Figure  27.  The  formula  above  gives 
sp=  0.7  s  and  ss=1.4  s.  These  estimates  are  comparable  to  the  standard  deviations  of  the 
residuals  for  Pn  and  Sn  arrival  times  shown  in  Figure  29  .  The  SDl’s  are  robust  estimates 
(include  about  98%  of  the  data)  and  the  SD2’s  include  all  data  i.e.,  outlying  observations. 

Some  of  the  scatter  in  the  arrival  time  residuals  above  may  be  a  result  of  the  uncertainty  in 
the  determinations  of  the  event  locations.  The  ground  truth  data  base  at  CSS  (Grant  and 
Coyne,  1992)  provide  epicenter  determinations  independent  of  the  IMS  locations.  Here 
we  use  observations  for  a  suite  of  earthquakes  near  Steigen  in  Northern  Norway  (see  map 
in  Figure  30 )  to  analyze  residuals  and  S/P  travel  time  ratios.  The  map  in  Figure  31  shows 
the  clustering  of  the  events  which,  with  the  exception  of  2  events,  are  within  10  km.  The 
epicenter  determinations  were  reported  by  Bergen  Seismological  Observatory  in  Norway 
and  were  based  on  measurements  at  a  local  network  of  5  stations. 

The  arrival  time  residuals  for  Pn  and  Sn  -  defined  as  observed  arrival  time  minus  that  com¬ 
puted  with  the  IMS  travel  time  model  and  the  origin  times  reported  by  Beigen  -shown  in 
Figure32  show  significant  scatter.  Apart  from  measurement  errors,  the  large  scatter  may  be 
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Figure  27:  Comparisons  of  origin  times  calculated  from  travel  time  ratios  of 
S  and  P  waves  and  the  associated  estimates  of  the  travel  time  ratios. 
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Figure  28:  Standard  deviation  of  origin  time  for  different  standard  deviations 
in  the  S/P.  from  0.1  to  1  per  cent,  as  a  function  of  epicentrai  distance.  The 
curves  were  obtained  with  a  Gaussian  random  number  generator  with  the 
standard  deviations  of  the  P  and  S  arrival  times  being  0.7  and  1.4  s,  respec¬ 
tively. 
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Figure  29:  Empirical  distribution  functions  of  Pn  and  Sn  arrival  time  residu¬ 
als  for  IMS  data.  The  lines  represent  approximations  to  Gaussian  distribu¬ 
tions;  the  dashed  and  solid  lines  correspond  to  robust  and  conventional 
estimates  of  the  standard  deviations  (SD1  and  SD2  respectively). 
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Figure  30:  Epicenters  of  the  events  that  were  analyzed  and  the  locations  of 
the  mini  arrays. 
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Figure  31:  Epicenters  of  the  Steigen  earthquakes  determined  by  the  Bergen 
Seismoiogical  Observatory. 
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Figure  32:  Comparison  of  Pn  and  Sn  residuals  of  the  Steigen  earth¬ 
quakes. 
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caused  by  errors  in  the  origin  times.  Except  for  two  outlying  data  points  each  for  ARAO 
and  NRAO,  the  Pn  and  Sn  residuals  for  a  given  array  also  correlate  as  one  might  expect 
with  the  derivative,  dT(Sn)/dT(Pn)=1.75,  indicated  as  a  dashed  line  in  the  diagram.  How¬ 
ever,  the  median  residuals  for  Pn  and  Sn  also  suggest  lateral  variations  in  Pn  and  Sn  as 
indicated  by  Table  4.  The  fastest  Pn  travel  time,  to  ARAO,  is  about  1.5s  less  than  that  to 
NRAO.  We  also  note  that  the  Sn  travel  times  to  these  arrays  are  about  2  s  slower  than  one 
would  expect  from  a  1.75  S/P  travel  time  ratio,  whereas  the  travel  time  of  Sn  to  FIAO  rela¬ 
tive  to  P  is  more  in  accordance  with  a  1.75  S/P  ratio. 

The  scatter  in  the  Pn  and  Sn  residual  is  also  reflected  in  the  S/P  travel  time  ratios,  which 
are  all  slightly  above  that  for  the  IMS  travel  time  tables  (Figure  33  ).  We  note  that  both 
NRAO  and  ARAO  appear  to  have  larger  ratios  than  does  FIAO. 

Examples  of  Depth  Estimation  Using  IMS  Data 

We  applied  the  LOCEX  program  to  events  shown  in  the  map  of  Figure  30,  both  with  a  free 
depth  and  with  the  origin  time  constrained  by  the  S/P  travel  time  ratio.  The  events  and  the 
resulting  depth  estimates  are  listed  in  Table  4.  The  events  include  the  suite  of  14  earth¬ 
quakes  near  Steigen  analyzed  above.  In  addition  there  are  two  groups  of  events,  the  nature 
of  the  source  of  which  has  not  been  firmly  established,  but  which  are  likely  to  be  earth¬ 
quakes  and  explosions  -  mining  or  under  water.  These  two  event  types  are  indicated  as 
possible  earthquakes  and  explosions  respectively  in  Table  4  and  are  marked  with  different 
symbols  in  Figure  30. 

The  point  estimates  of  the  two  types  of  depth  estimates  are  compared  in  Figures  34  and  35. 
The  scatter  diagram  in  Figure  34  shows  that  the  two  types  of  estimate  generally  agree,  but 
there  are  also  differences.  Several  of  the  events  designated  possible  explosions  are  placed 
at  depths  below  1  km  by  the  free  depth  method.  Furthermore,  there  are  several  possible 
earthquakes  that  have  point  estimates  at  zero  depth  for  the  free  depth  method.  Apparent 
inconsistencies  like  these  can  be  noted  for  depth  estimates  based  on  constrained  origin 
time  as  well,  but  they  appear  less  frequently  as  demonstrated  by  the  histograms  in  Figure 
35,  where  the  empirical  distributions  of  the  depth  estimates  of  the  two  kinds  of  events  and 
the  two  methods  are  compared. 

The  differences  observed  in  the  two  kinds  of  depth  estimates  is  no  evidence  that  depth 
from  constrained  origin  time  would  generally  provide  the  more  accurate  depth  for  several 
reasons.  For  one  thing,  the  data  analyzed  here  is  limited  to  a  fairly  small  number  of  events, 
the  exact  depths  of  which  are  not  known.  Furthermore,  the  uncertainties  of  the  point  esti¬ 
mates  have  not  been  analyzed  here  as  that  would  require  some  consideration  beyond  stan¬ 
dard  statistical  practice  for  the  depths  based  on  constrained  origin  times.  Nevertheless  it  is 
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Figure  33:  Empirical  distribution  functions  of  the  S/P  travel  time  ratios  to 
ARAO,  FIAO,  and  NRAO  for  the  Steigen  earthquakes. 
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Figure  34:  Comparison  of  depth  estimates  obtained  with  a  free  depth  and 
with  origin  time  constrained  with  the  S/P  travel  time  ratio. 
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Figure  35:  Histograms  of  depths  for  possible  earthquakes  and  mining 
events  determined  with  conventional  free  depth  method  and  with  origin 
times  constrained  with  S/P  travel  time  ratio. 
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fell  dial  the  results  warrant  further  analyses  and  comparisons  of  the  qualities  of  depth  esti¬ 
mates  based  on  constrained  origin  times. 

Concluding  Remarks 

In  this  note  we  have  used  simple  simulations  and  calculations  on  actual  data  to  compare 
depth  estimates  obtained  by  standard  location  algorithms  with  those  based  on  origin  times 
constrained  with  S/P  travel  time  ratio.  Because  of  die  limited  scope  of  the  simulations  and 
the  data  analyzed  the  observed  improvement  in  the  depth  estimates  based  on  constrained 
origin  times  has  to  be  confirmed  with  more  comprehensive  evaluations.  The  variability  in 
the  estimates  of  S/P  ratio  for  Fennoscandia  also  suggest  that  to  be  truly  useful,  the  regional 
dependence  of  has  to  be  mapped. 
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Table  4:  Event  source  parameters  and  depth  estimates 


IMS  Hypoc enter  Solution 

Depths  (km) 

Type  of 
event 

Date 

origin  time 

LaL 

Loo. 

Depth 

ML 

Free 

Const 
r.  OT 

90/01/06 

22:55:42 

6318 

22.67 

0.0 

132 

3.9 

3.8 

Foss.  eq. 

90/03/29 

04:11:19 

62.07 

6.06 

0.0 

2.30 

12 

4.7 

Foss.  eq. 

90/05/24 

56.12 

12.82 

0.0 

132 

26.9 

21.1 

Foss.  eq. 

90/08/22 

04:08:30 

63.86 

20.81 

0.0 

2.12 

11.0 

10.6 

Foss.  eq. 

90/10/28 

02:26:10 

66.85 

2430 

0.0 

1.97 

18.8 

8.8 

Foss.  eq. 

90/10/30 

23:08:11 

67.65 

15.41 

0.0 

2.29 

91 

10.0 

Foss.  eq. 

91/04/17 

17:47:29 

67.11 

24.16 

0.0 

2.11 

0.0 

41 

Foss.  eq. 

91/06/13 

10:48:23 

67.75 

20.16 

0.0 

149 

0.0 

0.0 

Foss.  eq. 

91/07/13 

01:42:21 

67.99 

23.86 

0.0 

119 

1.7 

0.3 

Foss.  eq. 

91/09/23 

64.52 

2134 

0.0 

3.15 

0.0 

El 

Foss.  eq. 

91/10/28 

13:14:34 

59S7 

6.28 

0.0 

175 

0.0 

16.3 

Foss.  eq. 

91/10/28 

16:21:27 

63.83 

20.78 

0.0 

115 

0.0 

0.0 

Foss.  eq. 

91/11/08 

01:20:36 

58.01 

14.84 

0.0 

121 

6.2 

0.0 

Foss.  eq. 

91/12/11 

16.-00:23 

66.95 

2436 

0.0 

119 

4.4 

4.3 

FOss.  eq. 

91/12/15 

15:18:48 

62.03 

17.99 

0.0 

174 

5.3 

5.0 

Foss.  eq. 

92/01/01 

08:03:58 

67.75 

14.84 

12.1 

4.50 

0.1 

2.9 

Steigen  eq. 

92/01/01 

08:39:02 

67.72 

14.85 

12.1 

1.60 

2.9 

1.5 

Stcigen  eq. 

92/01/04 

03:43:43 

67.71 

14.93 

0.0 

2.00 

15 

0.0 

Steigen  eq. 

92/01/04 

04:15.-04 

67.71 

14.90 

12.1 

3.60 

5.6 

6.9 

Steigen  eq. 

92/01/04 

06:00:53 

67.70 

14.87 

12.1 

wxm 

7.5 

5.4 

Steigen  eq. 

92/01/04 

09:06:31 

67.75 

14.92 

12.1 

150 

4.4 

5.3 

Steigen  eq. 

92/01/05 

01:20:50 

67.73 

14.78 

15.5 

1.70 

5.2 

3.9 

Steigen  eq. 

92/01/05 

05:11:56 

67.71 

14.99 

0.0 

IQ! 

10.6 

12.4 

Steigen  eq. 
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Ifcble  4:  Event  aource  parameters  and  depth  estimates 


IMS  Hypocentcr  Solution 

Depths  (km) 

Type  of 
event 

Date 

origin  time 

Lat 

Loo. 

Depth 

ML 

Free 

Const 
r.  OT 

92/01/06 

08:26:22 

67.70 

14.96 

gETil 

4.6 

32 

Steigen  eq. 

92/01/10 

07:46:36 

6727 

23.63 

0.0 

244 

73 

132 

Foss.  eq. 

92/01/10 

22:2S:42 

67.71 

14.98 

22 

-9.99 

28 

22 

Steigen  eq. 

92/01/25 

11:57:34 

67.69 

14.93 

0.0 

0.4 

0.0 

Steigen  eq. 

92/01/25 

12:16:48 

67.71 

14.88 

121 

5.9 

28 

Steigen  eq. 

92/01/25 

12:26:29 

67.73 

15.11 

0.0 

62 

3.0 

Steigen  eq. 

92/01/25 

19:12:52 

67.74 

14.61 

121 

43 

3.8 

Steigen  eq. 

92/02/19 

06:39:33 

59.24 

10.89 

24.5 

3.26 

28.3 

23.9 

Foss.  eq. 

92/04/16 

10:21:43 

64.23 

20.43 

0.0 

234 

0.0 

0.7 

Foss.  eq. 

92/06/10 

08:22:07 

58.98 

1833 

0.0 

238 

0.0 

0.0 

Poss.  expl. 

92/06/10 

11:46:33 

69.32 

3036 

0.0 

271 

0.0 

0.0 

Foss.  expl. 

92/06/19 

11:35:40 

6929 

3036 

0.0 

266 

0.0 

0.0 

Foss.  expl. 

92/06/20 

23:01:32 

6736 

15.68 

6.8 

274 

02 

0.8 

Foss.  eq. 

92/06/23 

12:21:33 

58.19 

18.96 

0.0 

250 

0.0 

0.0 

Poss.  expl. 

92/06/23 

12:36:11 

5732 

19.20 

0.0 

260 

0.0 

0.0 

Poss.  expl 

92/06/23 

12:40:30 

57.73 

18.99 

0.0 

2.51 

0.0 

0.0 

Poss.  expl. 

92/06/28 

12:23:33 

61.74 

5.29 

0.0 

2.62 

23 

4.9 

Foss.  eq. 

92/07/01 

07:12:23 

59.88 

30.45 

0.0 

254 

0.0 

0.0 

Poss.  expl. 

92/07/03 

11:32:54 

69.37 

30.71 

0.0 

252 

0.8 

0.0 

Poss.  expl. 

92/07/11 

14:32:28 

62.13 

4.03 

7.3 

264 

15.0 

27.4 

Foss.  eq. 

92/08/03 

11:16:02 

66.32 

13.58 

11.8 

2.69 

11.8 

10.7 

Foss.  eq. 

92/08/06 

07:32:41 

59.85 

6.30 

4.7 

3.15 

3.7 

32 

Foss.  eq. 

92/08/28 

18:37:52 

59.94 

6.62 

8.6 

2.63 

2.9 

3.2 

Foss.  eq. 

I 
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Table  4;  Event  source  parameters  and  depth  estimates  _ 


Hypocentcr  Solution 

Depths  (km) 

TVpe  of 

Lat. 

Lon. 

Depth 

ML 

Free 

Const 
r.  OT 

event 

92/09/23 

11:37.-02 

6938 

30.73 

0.0 

162 

1.7 

0.0 

Pou.  expL 

12331 

12:22:51 

6930 

30.77 

0.0 

2.86 

0.0 

Foss.  expL 

92/10/17 

02:19:38 

62.14 

2.65 

0.0 

198 

16 

3.6 

Foss.  eq. 

92/10/23 

14:50:23 

6938 

34.92 

0.0 

154 

3  2 

0.0 

Foss.  expL 

92/10/30 

12:5433 

6939 

30.73 

0.0 

178 

0.0 

0.0 

Foss.  expL 

92/11/18 

10:20:08 

64.61 

31.07 

0.0 

169 

0.0 

0.0 

Poss.  expL 

92/11/20 

69.29 

30.61 

0.0 

189 

6.4 

0.0 

Pass.  expL 

92/11/22 

04:46:45 

67.61 

33.47 

0.0 

157 

0.0 

0.0 

Foss.  expL 

92/11/30 

14:10:39 

69.42 

31.64 

0.0 

190 

0.0 

0.0 

Foss.  expL 

92/12/11 

13:16:31 

6931 

30.62 

0.0 

159 

0.0 

0.0 

Poss.  expL 

92/12/18 

12:09:52 

68.11 

33.13 

0.0 

162 

33 

0.0 

Poss.  expL 

92/12/23 

12:29:41 

69.43 

30.79 

0.0 

169 

0.0 

0.0 

Poss.  expL 

92/12/30 

11:30:59 

67.94 

3230 

0.0 

155 

113 

1.4 

Poss.  expL 

92/12/30 

15.-05:25 

69.34 

30.65 

0.0 

174 

0.0 

0.0 

Poss.  expL 

93/01/05 

10:19:35 

64.71 

17.17 

1.8 

3.66 

0.0 

3.8 

Foss.  eq. 

21:44:08 

73.17 

6.30 

0.0 

182 

1.7 

1.7 

Foss.  eq. 

93/01/11 

04:11:48 

57.95 

7.58 

0.0 

175 

0.4 

03 

Poss.  expL 

93/01/15 

12:28:23 

69.35 

30.63 

0.0 

150 

0.0 

0.0 

Poss.  expL 

93/01/16 

12:30:11 

68.12 

33.07 

0.0 

167 

1.4 

0.0 

Poss.  expL 

93/01/29 

12:40:20 

69.38 

30.69 

0.0 

2.82 

0.0 

0.0 

Poss.  expL 

93/01/30 

13:07:09 

67.62 

33.76 

0.0 

2.53 

0.0 

0.0 

Poss.  expL 

93/02/05 

12:13:38 

69.39 

30.74 

0.0 

164 

0.1 

0.0 

Poss.  expL 

93/02/05 

16:07:56 

70.71 

14.92 

0.0 

2.58 

0.0 

0.0 

Poss.  eq. 

Assessing  Depth  Estimation  Capability:  Examples  for 
Global  Networks 

Hans  Israelsson 


Introduction 

Estimating  the  depth  of  seismic  events  is  one  important  function  of  station  networks. 
Depth  is  an  important  discriminant  for  differentiating  between  explosions  and  earthquakes 
and  it  is  an  important  parameter  in  studies  of  seismic  activity  and  associated  geophysical 
processes.  Yet  the  capability  to  estimate  depths  of  seismic  networks  are  rarely  established 
and  station  networks  are  rarely  designed  to  allow  depth  estimation  in  an  optimum  manner. 

In  this  note  we  perform  some  computational  experiments  to  outline  the  depth  estimation 
capability  of  the  GSETT-2  global  network  and  a  hypothetical  network  consisting  primarily 
of  high  sensitive  array  stations.  We  also  discuss  the  use  of  stations  in  the  oceanic  areas  to 
supplement  the  hypothetical  network  to  enhance  its  capabilities  of  depth  estimation  in  the 
oceans. 

The  Networks 

Simulations  of  the  detection  capability  of  the  GSETT-2  network  (see  Bratt,  1992)  have 
been  carried  out  by  Bratt  (1993)  and  are  summarized  in  the  map  of  Figure  36a.  The  con¬ 
tours  were  obtained  by  the  NetSim  software  package  (Sereno  et  al.,  1990)  and  parameters 
developed  and  employed  by  Bratt  (1993).  The  capabilities  can  be  compared  with  that  of  a 
hypothetical  network  of  SO  sensitive  stations  (ALPHA  network)  in  Figure  36b  (also  stud¬ 
ied  by  Bratt,  1993).  Finally,  in  Figure  36c  are  the  contours  of  the  detection  threshold  for 
another  hypothetical  network  consisting  of  the  ALPHA  network  supplemented  with  50 
oceanic  stations,  selected  from  sites  of  existing  or  previously  operated  seismic  stations  in 
oceanic  areas,  with  noise  and  detection  characteristics  equal  to  those  of  a  typical  oceanic 
station  (SNZO)  of  the  GSETT-2  network.  The  potential  of  oceanic  stations  for  explosion 
monitoring  have  earlier  been  analyzed  and  discussed  by  Phillips  and  McCowan  (1978). 
Median  thresholds,  ranges,  and  standard  deviation  are  given  for  different  geological 
regions  -  oceanic,  tectonic,  shield/platform  -  following  the  regionalization  by  Jordan 
(1981).  We  are  systematically  using  medians  instead  of  mean  values  since  grid  points  for 
which  capabilities  are  calculated  may  be  close  to  some  stations  and  therefore  will  be  low. 
Similarly,  we  use  a  robust  estimate  of  the  standard  deviation  that  is  insensitive  to  outlying 
data  values.  The  ALPHA  network  has  a  more  uniform  detection  capability  although  the 
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Figure  36:  The  90%  detection  thresholds  of  the  GSETT-2  network  (a),  the 
ALPHA  network  (b),  and  a  network  consisting  of  the  ALPHA  staitons  and  50  oce¬ 
anic  stations  (c). 
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detection  threshold  is  about  two  tenths  of  a  unit  higher  in  oceanic  ar  “as  than  in  other  areas. 

Station  detection  thresholds  for  the  teleseismic  detection  window  at  the  50%  level  are 
shown  in  the  histograms  of  Figure  37.  They  were  obtained  as  the  median  thresholds  from 
NetSim  calculations  with  one  station  networks.  Since  the  medians  represent  detection 
thresholds  at  90  degrees  distance,  a  factor  of  0.35  magnitude  units  (corresponding  to  the 
attenuation  difference  between  90  and  60  degrees,  60  degrees  being  in  the  center  of  the 
teleseismic  detection  window)  were  subtracted.  These  thresholds  showed  close  agreement 
with  thresholds  for  the  GSETT-2  station  distribution  obtained  from  comparisons  of 
detected  and  undetected  events  reported  by  the  NEIS.  The  mb(NEIS)  were  corrected  for 
bias  using  empirical  corrections  from  the  results  of  Ringdal  (1986).  The  distribution  of  the 
detection  thresholds  for  the  ALPHA +OCEANIC  network  is  dominated  by  the  50  oceanic 
stations.  Phillips  and  McGowan  (1979)  refer  to  observations  of  low  noise  levels  in  one 
early  ocean  seismometer  experiment.  Subsequent  studies  have,  however,  not  substantiated 
these  early  observations  (Hedlin  and  Orcutt,  1989).  Hence  the  large  noise  amplitude  and 
comparatively  high  detection  threshold  for  the  oceanic  stations. 

Depth  Estimation  and  Simulations 

There  are  several  techniques  for  estimating  focal  depth,  the  most  commonly  employed  of 
which  is  the  estimation  resulting  from  U»e  standard  location  of  first  arrival  times  of  P 
phases.  It  is  common  knowledge  that  such  estimates  are  usually  associated  with  large 
uncertainties  unless  the  variance  of  the  travel  time  derivatives  with  respect  to  depth  of  the 
recording  stations  is  large.  Depth  estimates  based  on  secondary  phases  often  can  often  dra¬ 
matically  improve  the  quality  of  the  estimate,  but  secondary  phases  are  not  always  avail¬ 
able. 

Depth  estimates  used  for  the  purpose  of  explosion  monitoring  do  not  necessarily  have  the 
same  requirements  as  those  for  precise  focal  depth  determinations  of  earthquakes.  For 
example,  assessing  depth  accuracy  based  on  S  wave  arrival  times  (one  of  the  most  effi¬ 
cient  depth  estimation  methods)  Gombeig  et  al.  (1990)  concluded  that  in  order  to  achieve 
an  accuracy  of  about  1.5  km  would  require  S  waves  recorded  within  1.4  times  the  focal 
depth.  Applying  this  criterion  for  world-wide  seismicity  would  require  a  network  with  a 
v  -ry  large  number  of  stations.  For  explosion  monitoring  it  is  merely  sufficient  to  establish 
that  the  focal  depth  is  deeper  than  some  limit,  say  10  km.  Here,  however,  we  will  look  at 
actual  depth  estimates  rather  than  upper  limits.  We  analyze  depth  estimation  based  jn  two 
kinds  of  secondary  phases  (1)  surface  reflected  pP  and  sP  recorded  at  teleseismic  distances 
and  (2)  S  waves  recorded  primarily  at  regional  distances. 
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Figure  37:  The  50%  detection  thresholds  of  the  teleseismic  window  for  station  in 
the  GSETT-2  and  ALPHA  networks  and  a  network  consisting  of  the  ALPHA  sta¬ 
tions  and  50  oceanic  stations. 
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Surface  Reflected  Phases  -  Moveout  and  Azimuthal  Coverage 

From  past  experience  it  appears  that  the  observability  of  depth  phases,  pP  and  sP,  from 
earthquakes  is,  in  many  instances  unpredictable,  and  is  therefore  difficult  to  model.  Simple 
models  for  the  radiation  pattern  of  double  couple  sources  of  short  duration,  1  s  ex-  so,  have 
in  several  cases  been  demonstrated  to  accurately  describe  the  excitation  and  recording  of 
such  these  phases.  It  appears,  however,  that  the  success  rate  of  this  model  has  not  been 
firmly  established  and  it  seems  that  in  a  large  number  of  cases,  other  factors  dominate  the 
character  of  recorded  P  seismograms  (Douglas  et  al.,  1990).  Therefore  we  limited  the 
modelling  of  the  accuracy  of  depths  determined  from  depth  phases  to  simple  factors  that 
from  experience  control  their  observability  to  some  degree  azimuthal  coverage  of  record¬ 
ing  stations,  and  the  degree  of  moveout  of  the  depth  phases  that  recording  stations  may 
provide.  Both  these  factors  are  determined  by  the  geometry  of  the  station  network.  Only 
stations  that  detect  the  event  in  the  teleseismic  window  (30-90  degrees)  are  assumed  to 
contribute  to  moveout  and  azimuthal  coverage.  This  is  somewhat  conservative  as  depth 
phases  may  be  detected  and  positively  identified  at  shorter  distances  in  spite  of  the  com¬ 
plexity  of  regional  wave  propagation. 

The  moveout  is  given  as  the  percentage  observed  of  the  total  moveout  possible  in  the  dis¬ 
tance  range  30  to  90  degrees  and  was  obtained  by  Monte  Carlo  simulations  similar  to 
those  employed  in  the  location  module  of  the  NetSim  program  (Sereno  et  al.,  1990).  For  a 
given  network  -  GSETT  and  ALPHA  -  the  moveout  provided  for  an  event  with  given  loca¬ 
tion  and  magnitude  represents  the  average  of  a  large  number  of  randomizations  of  signal 
and  noise  amplitudes  assuming  the  50%  detection  thresholds  shown  in  Figure  37.  The 
results  are  summarized  in  Figure  38  for  the  GSETT  and  the  ALPHA  networks;  the 
ALPHA  and  ALPHA +OCEAN  network  provided  very  similar  results.  For  each  of  the  two 
networks  in  Figure  38,  the  median  of  the  percent  moveout  is  shown  as  a  function  of  mag¬ 
nitude.  Medians  are  given  according  to  type  of  geological  region  -  oceanic,  tectonic,  and 
shield/platform.  The  median  of  the  detection  thresholds  for  ?  three  types  of  geological 
region  (see  Figure  36a-b)  are  also  indicated  on  the  diagrams  as  a  shaded  zone  for  the  sake 
of  comparison.  We  note  that  the  ALPHA  network  reaches  close  to  100%  moveout  at  about 
magnitude  4.0  and  only  at  the  low  end  of  the  range  for  the  detection  thresholds  does  the 
moveout  go  below  50%.  The  trend  for  the  GSETT-2  network  relative  to  its  detection 
thresholds  is  rather  similar,  except  that  even  at  larger  magnitudes,  100%  moveout  is  not 
attained. 

In  Figure  39  we  compare  similar  curves  for  the  azimuthal  coverage  given  as  the  percent¬ 
age  of  8, 45  degree  azimuth  sectors,  in  which  the  event  was  detected.  The  results  in  Figure 
39  were  obtained  with  simulations  similar  to  those  used  for  the  moveout.  For  each  simula¬ 
tion,  the  azimuthal  coverage  of  the  detecting  stations  was  calculated  by  counting  the  num- 


81 


GSETT-2 


MEDIAN  OF  MOVEOUT 


ALPHA 


Figure  38:  Median  of  the  moveout  observed  of  the  total  possible  moveout  in  the 
distance  range  30  to  90  degrees  (in  per  cent)  as  a  function  of  ganitude,  mb,  or 
median  of  the  detection  thresholds  (see  Fig  36a-b) 
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the  ALPHA  network.  The  shaded  zones  represent  the  median  of  the  detection 
thresholds  (see  Fig  36a-b). 


ber  of  45  degree  azimuthal  sectors  around  the  source  to  station  azimuth  that  woe  covered 
by  one  or  more  detecting  stations.  Centering  of  the  azimuth  at  the  median  of  the  station 
azimuth  was  done  prim-  to  the  counting  in  order  to  avoid  stations  clustered  in  one  45 
degree  sector  being  counted  as  two  sectors.  The  50%  coverage  is  indicated  as  a  dashed 
line  in  the  diagrams.  Analysis  of  G  SETT-2  indicated  that  50%  azimuthal  coverage 
improves  the  probability  of  detecting  and  positively  identifying  depth  phases.  The  curves 
in  Figure  39  for  the  GSETT-2  network  do  not  reach  this  level  until  well  over  half  a  magni¬ 
tude  unit  above  its  detection  thresholds,  whereas  this  gap  for  the  ALPHA  network  is  less 
pronounced,  but  still  significant  This  means  that  for  events  between  the  detection  thresh¬ 
old  and  the  50%  level,  the  probability  of  detecting  depth  phases  will  be  reduced. 

Depth  Errors  Using  S  Phases 

The  use  of  S  wave  arrival  times  for  the  depth  estimation  is  well  established  and  its  short¬ 
comings  and  capabilities  lend  themselves  more  easily  to  modelling  than  do  surface 
reflected  phases.  We  employed  the  NetSim  program  with  simulation  parameters  used  for 
GSE  by  Bratt  (1993)  and  describe  the  depth  error  with  90%  confidence  limits  obtained  by 
Monte  Carlo  simulations.  In  order  to  illustrate  the  effect  of  S  waves  on  the  estimated  depth 
error  we  compare  the  ratio  of  the  depth  errors  from  determinations  with  only  P  waves  and 
with  both  P  and  S  waves  in  Figure  40.  The  event  magnitude  was  set  to  mb=4.0  and  the 
depth  to  25  km.  There  is,  as  one  should  expect,  a  general  improvement  -  the  ratio  is  larger 
than  1.0.  However,  there  is  also  significant  variation  in  the  ratio;  the  median  for  the  differ¬ 
ent  geological  regions  varies  from  about  20  to  100%,  but,  as  indicated  by  the  ranges, 
improvements  can  be  over  500%. 

Figures  41  and  42  summarize  simulations  of  the  depth  error  for  depth  estimates  based  on  S 
waves  as  a  function  of  mb  for  the  three  networks.  The  event  is  assumed  to  be  at  25  km 
(indicated  with  a  dashed  horizontal  line)  and  curves  for  the  median  and  upper  limit  are 
shown  in  Figures  41  and  42,  respectively.  As  in  Figure  38  and  39,  the  range  of  the  median 
detection  thresholds  (from  Figure  36)  are  included  in  the  diagrams  for  the  sake  of  compar¬ 
ison.  The  curves  for  the  ALPHA  and  the  ALPHA+OCEAN  networks  in  Figure  41  and  42 
are  quite  similar  suggesting  very  little  effect  of  the  oceanic  stations  on  the  depth  error  for 
events  with  magnitudes  below  5.0.  In  order  to  achieve  improvement  of  the  ALPHA+- 
OCEAN  network  relative  to  the  ALPHA  network,  the  detectability  of  the  oceanic  stations 
have  to  be  enhanced. 

For  all  three  networks  there  is  a  clear  gap  between  the  depth  estimation  capability  as 
described  by  the  90%  error  limits  and  the  detection  thresholds.  We  define  the  gap  as  the 
difference  between  the  magnitude  at  which  the  90%  error  equals  25  km  (the  actual  depth 
of  the  event)  and  the  90%  detection  threshold.  The  gap  is  most  pronounced  for  oceanic 
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Figure  40:  The  ratio  of  depth  errors  from  determinations  with  only  P  waves  and 
with  both  P  and  S  waves. 


Figure  41:  Median  of  depth  error  as  a  function  of  mb  for  estimates  based  on  P 
and  S  waves  for  three  networks:  GSETT-2,  ALPHA,  and  ALPHA  +  a  50  station 
oceanic  network.  The  shaded  zones  represent  the  median  of  the  detection  thresh¬ 
olds  (see  Fig  36a-b). 


Figure  42:  Upper  limit  of  depth  error  as  a  function  of  mb  for  estimates  based  on  P 
and  S  waves  for  three  netowrks:  GSETT-2,  ALPHA,  and  ALPHA  +  a  50  station 
oceanic  network.  The  shaded  zones  represent  the  median  of  the  detection  thresh¬ 
olds  (see  Fig  36a-b). 


structures  for  the  GSETT  network  where  it  is  almost  one  magnitude  unit  For  other  cases  it 
stays  around  0.5  magnitude.  The  results  in  Figure  42,  showing  the  upper  limits  of  the 
depth  error,  show  that  the  differences  between  oceanic  and  other  structures  are  less  pro¬ 
nounced  than  they  are  for  the  median  value  results  in  Figure  41.  The  gaps  have,  of  course, 
widened  compared  to  the  median  errors.  For  the  GSETT  network,  the  increase  in  the  gap 
is  dramatic,  whereas  for  the  ALPHA  network  it  increases  about  0.3  magnitude  units  for 
tectonic  and  shield/platform  structure.  The  upper  limit  is  perhaps  conservative  as  a  mea¬ 
sure  of  the  capability,  but  a  large  upper  limit  also  reflects  that  there  are  areas  with  poor 
capability  even  if  the  median  levels  suggest  a  high  capability. 

One  reason  for  the  gap  between  detection  threshold  and  depth  estimation  capability  may 
be  due  to  the  apparent  lower  detection  capability  of  S  waves  as  illustrated  by  the  example 
in  Figure  43.  Here  the  difference  in  detection  thresholds  (90%)  for  detecting  both  P  and  S 
waves,  and  only  P  waves  at  three  stations  is  contoured.  As  S  waves  are  often  difficult  to 
identify  and  properly  time  the  requirement  of  detecting  S  waves  at  three  stations  rather 
than  just  the  minimum  one  was  used.  The  median  difference  in  thresholds  in  Figure  43  is  a 
few  tenths  of  a  magnitude  unit,  but,  as  indicated  by  the  ranges,  can  reach  higher  values 
close  to  a  magnitude  unit. 

The  gap  between  detection  and  depth  estimation  capabilities  means  that  some  events  will 
have  depths  that  are  only  poorly  determined.  In  order  to  illustrate  this  effect  we  applied  the 
depth  estimation  simulations  to  the  ALPHA  network,  not  to  events  as  a  function  of  magni¬ 
tude  on  an  equally  spaced  grid,  as  above,  but  to  events  that  model  the  world  wide  seismic 
activity.  For  this  purpose  we  generated  a  sample  of  earthquakes  corresponding  to  5  days  of 
activity  with  the  simulation  program  by  North  (1983).  This  sample  included  2997  events 
with  magnitudes  down  to  mb=2.58.  For  each  event,  we  then  calculated  the  depth  error 
using  a  large  number  of  Monte  Carlo  simulations.  The  resulting  depth  errors  are  compared 
in  Figure  44  with  the  actual  depths  of  the  synthetic  events.  Only  events  that  were  detected 
with  90%  probability  are  included  in  the  diagrams  of  Figure  44.  In  the  left  and  right 
frames  we  show  the  differences  between  the  90%  error  and  depth  vs.  magnitude  and  actual 
depth,  respectively.  Separate  symbols  -  plus  signs  and  open  circles  -  are  used  for  events 
with  the  difference  less  than  and  greater  than  0  indicated  as  a  dashed  line.  This  separation 
maps  events  with  acceptable  and  poorly  determined  depths  relative  to  the  actual  depths. 
Events  below  the  line  would  be  “screened”  as  earthquakes,  whereas  data  points  above  the 
dashed  line  are  too  large  to  allow  such  a  classification.  A  screening  level  of  0  km  may  be 
too  low  for  monitoring  purposes  and  it  may  be  more  prudent  to  set  it  at,  say  10  km.  We 
note  from  the  scatter  diagram  to  the  right  that  events  with  actual  depths  shallower  than  10 
km  are  generally  not  “screened”  as  earthquakes.  This  is  desirable  as  they  represent  a  depth 
shallow  enough  for  at  least  theoretically  be  explosions.  The  10  km  depth  line  is  indicated 
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Figure  43:  The  difference  in  detection  thresholds  (90%)  for  detecting  both  P  and  S 
waves,  and  only  P  waves  at  three  stations. 
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Figure  43:  The  difference  in  detection  thresholds  (90%)  for  detecting  both  P  and  S 
waves,  and  only  P  waves  at  three  stations. 


in  the  right  frame,  which  also  gives  the  number  of  events  at  depths  below  10  km  that  are 
not  screened,  135  or  27  events/day.  This  corresponds  to  about  25%  of  the  total  number  of 
events  that  were  detected  with  a  90%  probability.  Although  the  actual  depth  strongly 
determines  this  number,  the  left  frame  also  indicates  that  these  are  primarily  small  events 
most  of  which  have  magnitudes  less  than  mb=4.0.  The  mb  detection  thresholds  (from  Fig¬ 
ure  36)  are  indicated  in  the  left  frame.  Finally,  the  geographical  distribution  of  the  events 
are  shown  in  Figure  45,  where  events  that  are  above  the  “screening”  line  and  at  depths 
below  10  km  are  marked  filled  circles.  There  is  no  clear  concentration  of  these  events. 

Concluding  Remarks 

The  simulations  and  examples  presented  here  suggest  that  the  ALPHA  network  has  not 
only  a  much  lower  and  more  geographically  uniform  detection  threshold  than  the  GSETT- 
2  network,  but  has  a  much  improved  depth  estimation  capability  as  well. 

The  simulations  also  suggest  that  only  marginal  improvements  would  be  obtained  by  add¬ 
ing  50  oceanic  stations  to  the  ALPHA  network,  unless  the  sensitivity  of  these  ocean  based 
stations  are  significantly  improved  as  compared  to  typical  levels  of  stations  in  operation, 
currently  and  in  the  past 

Although  the  observability  of  depth  phases  is,  in  large  degree,  unpredictable  and  difficult 
to  model  generally,  it  appears  that  the  ALPHA  network  is  well  distributed  to  observe  such 
phases  in  the  teleseismic  window  with  regard  to  moveout  of  depth  phases,  as  well  as  to 
azimuthal  distribution  around  event  azimuths.  These  qualities  are  reduced,  however,  near 
the  detection  threshold,  and  there  is  a  gap  of  a  few  tenths  of  a  magnitude  unit  where  the 
probability  of  observing  depth  phases  appears  low. 

There  is  a  similar  gap  for  depth  estimates  based  on  P  and  S  wave  arrival  times.  On  the 
average,  this  gap  does  not  appear  as  pronounced  for  events  in  tectonic  and  shield/platform 
areas.  The  simulations  indicate  that  the  gap  is  partly  caused  by  the  lower  detection  capa¬ 
bility  of  S  waves  of  the  ALPHA  network  and  therefore  may  be  narrowed  with  enhanced 
detectability  of  S  waves.  The  significance  of  narrowing  the  gap  is  illustrated  with  exam¬ 
ples  of  synthetically  generated  events.  A  fairly  large  percentage  of  a  large  daily  average  of 
detected  and  located  events  at  depths  below  10  km  would  otherwise  have  poorly  con¬ 
strained  depths.  One  important  function  of  the  BETA  network  of  the  GSE  system  would 
be  to  strengthen  observations  of  S  waves  for  depth  determinations.  The  GAMMA  network 
could  also  contribute  in  this  regard,  but  its  role  in  the  GSE  system  is  currently  not  clearly 
defined. 
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